!/---------------------------------------------------------------------------/
! CVS VERSION INFORMATION
! $Id$
! $Name$
! $Revision$
!/===========================================================================/

MODULE MOD_RRKF_OBS
# if defined (RRKF)
  USE RRKVAL
  USE MOD_RRK
  USE CONTROL
  USE MOD_NCTOOLS
  USE MOD_UTILS
  USE MOD_PREC
  USE MOD_WD
  USE ALL_VARS, only: vxmin,vymin
  USE MOD_STARTUP
  USE MOD_INPUT

  IMPLICIT NONE
  SAVE

  REAL(SP),ALLOCATABLE :: AW0G(:,:),AWXG(:,:),AWYG(:,:)

!--Data Assimilation Parameters for Current Nudging  Assimilation
!

  INTEGER  :: CUR_MAX_LAYER          !!MAXIMUM NUMBER OF VERTICAL DATA FROM ANY OBS POINT


!
  !
!
!
!--Current Assimilation Object Type
!
  TYPE ASSIM_OBJ_CUR
   INTEGER  :: N_TIMES                     !!NUMBER OF DATA TIMES
   INTEGER  :: N_INTPTS                    !!NUMBER OF INTERPOLATION POINTS 
   INTEGER  :: N_T_WEIGHT                  !!DATA TIME FOR CURRENT OBSERVATION WEIGHTING
   INTEGER  :: N_LAYERS                    !!NUMBER OF OBSERVATIONS IN THE VERTICAL
   INTEGER  :: N_CELL                      !!CELL NUMBER OF OBSERVATION
   REAL(SP) :: X,Y                         !!X AND Y COORDINATES OF OBSERVATION
   REAL(SP) :: T_WEIGHT                    !!TIME WEIGHT
   REAL(SP) :: DEPTH                       !!DEPTH AT OBSERVATION STATION (MOORING)
   REAL(SP) :: SITA                        !!LOCAL ISOBATH ANGLE AT OBSERVATION
   REAL(SP),ALLOCATABLE,DIMENSION(:)  :: ODEPTH      !!OBSERVATION DEPTHS
   INTEGER                        :: chose_flag      !!0 IS NO, 1 IS YES
   INTEGER                        :: INNODE     !! OBS NEAREST NODE
   REAL(DP),ALLOCATABLE,DIMENSION(:)  :: obs_days       !!DATA TIMES
   TYPE(TIME),ALLOCATABLE,DIMENSION(:)  :: TIMES       !!DATA TIMES
   REAL(SP),ALLOCATABLE,DIMENSION(:,:):: UO,VO       !!OBSERVATION DATA FOR X,Y VELOCITY COMPONENTS
   INTEGER,ALLOCATABLE,DIMENSION(:)   :: INTPTS      !!POINTS USED TO INTERPOLATE TO OBSERVATION LOC 
   INTEGER,ALLOCATABLE,DIMENSION(:,:) :: S_INT       !!SIGMA INTERVALS SURROUNDING CURRENT MEASUREMENT
   REAL(SP),ALLOCATABLE,DIMENSION(:,:):: S_WEIGHT    !!SIGMA WEIGHTING                               
   REAL(SP),ALLOCATABLE,DIMENSION(:)  :: X_WEIGHT    !!SPATIAL WEIGHTING FOR INTERPOLATION POINTS
  END TYPE ASSIM_OBJ_CUR

!
!--Temperature/Salinity Assimilation Object Type
!
   TYPE ASSIM_OBJ_TS  
     INTEGER  :: N_TIMES                                  !!NUMBER OF DATA TIMES
     INTEGER  :: N_INTPTS                                 !!NUMBER OF INTERPOLATION POINTS
     INTEGER  :: N_T_WEIGHT                               !!DATA TIME FOR CURRENT OBSERVATION WEIGHTING
     INTEGER  :: N_LAYERS                                 !!NUMBER OF OBSERVATIONS IN THE VERTICAL
     INTEGER  :: N_CELL                                   !!CELL NUMBER OF OBSERVATION
     REAL(SP) :: X,Y                                      !!X AND Y COORDINATES OF OBSERVATION
     REAL(SP) :: T_WEIGHT                                 !!TIME WEIGHT
     REAL(SP) :: DEPTH                                    !!DEPTH AT OBSERVATION STATION (MOORING)
     REAL(SP) :: SITA                                     !!LOCAL ISOBATH ANGLE AT OBSERVATION
   REAL(SP),ALLOCATABLE,DIMENSION(:)  :: ODEPTH      !!OBSERVATION DEPTHS
   INTEGER                        :: chose_flag      !!0 IS NO, 1 IS YES
   INTEGER                        :: INNODE      !! OBS NEAREST NODE
   REAL(DP),ALLOCATABLE,DIMENSION(:)  :: obs_days       !!DATA TIMES
     TYPE(TIME),ALLOCATABLE,DIMENSION(:)  :: TIMES        !!DATA TIMES
     REAL(SP), ALLOCATABLE, DIMENSION(:,:) :: TEMP        !!OBSERVATION DATA FOR TEMPERATURE 
     REAL(SP), ALLOCATABLE, DIMENSION(:,:) :: SAL         !!OBSERVATION DATA FOR SALINITY 
     INTEGER,  ALLOCATABLE, DIMENSION(:)   :: INTPTS      !!POINTS USED TO INTERPOLATE TO OBSERVATION LOC
     INTEGER,  ALLOCATABLE, DIMENSION(:,:) :: S_INT       !!SIGMA INTERVALS SURROUNDING CURRENT MEASUREMENT
     REAL(SP), ALLOCATABLE, DIMENSION(:,:) :: S_WEIGHT    !!SIGMA WEIGHTING
     REAL(SP), ALLOCATABLE, DIMENSION(:)   :: X_WEIGHT    !!SPATIAL WEIGHTING FOR INTERPOLATION POINTS
   END TYPE ASSIM_OBJ_TS    
!--EL ASSIMILATION Assimilation Object Type
!
   TYPE ASSIM_OBJ_EL
     INTEGER  :: N_TIMES                                  !!NUMBER OF DATA TIMES
     INTEGER  :: N_INTPTS                                 !!NUMBER OF INTERPOLATION POINTS
     INTEGER  :: N_T_WEIGHT                               !!DATA TIME FOR CURRENT OBSERVATION WEIGHTING
     INTEGER  :: N_LAYERS                                 !!NUMBER OF OBSERVATIONS IN THE VERTICAL
     INTEGER  :: N_CELL                                   !!CELL NUMBER OF OBSERVATION
     REAL(SP) :: X,Y                                      !!X AND Y COORDINATES OF OBSERVATION
     REAL(SP) :: T_WEIGHT                                 !!TIME WEIGHT
     REAL(SP) :: DEPTH                                    !!DEPTH AT OBSERVATION STATION (MOORING)
     REAL(SP) :: SITA                                     !!LOCAL ISOBATH ANGLE AT OBSERVATION
   REAL(SP),ALLOCATABLE,DIMENSION(:)  :: ODEPTH      !!OBSERVATION DEPTHS
   INTEGER                        :: chose_flag      !!0 IS NO, 1 IS YES
   INTEGER                        :: INNODE      !! OBS NEAREST NODE
   REAL(DP),ALLOCATABLE,DIMENSION(:)  :: obs_days       !!DATA TIMES
     TYPE(TIME),ALLOCATABLE,DIMENSION(:)  :: TIMES        !!DATA TIMES
     REAL(SP), ALLOCATABLE, DIMENSION(:,:) :: EL        !!OBSERVATION DATA FOR TEMPERATURE
     INTEGER,  ALLOCATABLE, DIMENSION(:)   :: INTPTS      !!POINTS USED TO INTERPOLATE TO OBSERVATION LOC
     INTEGER,  ALLOCATABLE, DIMENSION(:,:) :: S_INT       !!SIGMA INTERVALS SURROUNDING CURRENT MEASUREMENT
     REAL(SP), ALLOCATABLE, DIMENSION(:,:) :: S_WEIGHT    !!SIGMA WEIGHTING
     REAL(SP), ALLOCATABLE, DIMENSION(:)   :: X_WEIGHT    !!SPATIAL WEIGHTING FOR INTERPOLATION POINTS
   END TYPE ASSIM_OBJ_EL

!  el assimilation variables
   INTEGER                          :: N_ASSIM_EL !!NUMBER OF CURRENT OBSERVATIONS
   TYPE(ASSIM_OBJ_EL), ALLOCATABLE :: EL0_OBS(:)  !!CU ASSIMILATION DATA OBJECTS

                                                                                                                       
!
!--Current Data Assimilation Variables
!
   INTEGER                          :: N_ASSIM_CUR !!NUMBER OF CURRENT OBSERVATIONS 
   TYPE(ASSIM_OBJ_CUR), ALLOCATABLE :: CUR_OBS(:)  !!CURRENT ASSIMILATION DATA OBJECTS

!
!--Salinity/Temperature Data Assimilation Variables
!
   INTEGER                                N_ASSIM_T   !!NUMBER OF TEMPERATURE OBSERVATIONS
   TYPE(ASSIM_OBJ_TS), ALLOCATABLE     :: T0_OBS(:)    !!TEMP ASSIMILATION DATA OBJECTS


   INTEGER                                N_ASSIM_S   !!NUMBER OF TEMPERATURE OBSERVATIONS
   TYPE(ASSIM_OBJ_TS), ALLOCATABLE     :: S0_OBS(:)    !!TEMP ASSIMILATION DATA OBJECTS
CONTAINS                                                                                                                       
   subroutine deallocate_obs_data
   implicit none
   IF (EL_ASSIM) THEN
      DEALLOCATE(EL0_OBS)
   END IF
   IF (UV_ASSIM) THEN
      DEALLOCATE(CUR_OBS)
   END IF
   IF (T_ASSIM) THEN
       DEALLOCATE(T0_OBS)
   END IF
   IF  (S_ASSIM)  THEN
      DEALLOCATE(S0_OBS)
   END IF

   end subroutine deallocate_obs_data
   subroutine filter_obs_data
   IMPLICIT NONE
   INTEGER :: I,J,K
   IF (EL_ASSIM) THEN
      CALL FILTER_EL_OBS_DATA(INTTIME)
   END IF
   IF (UV_ASSIM) THEN
      CALL FILTER_UV_OBS_DATA(INTTIME)
   END IF
   IF (T_ASSIM) THEN
      CALL FILTER_T_OBS_DATA(INTTIME)
   END IF
   IF  (S_ASSIM)  THEN
      CALL FILTER_S_OBS_DATA(INTTIME)
   END IF 
   end subroutine filter_obs_data
   subroutine filter_el_obs_data(current_time)
   IMPLICIT NONE
   real(dp) :: time1
   TYPE(TIME) :: current_Time
!!-------------------------------------------------
   INTEGER I,J,K,NN,ECNT,ITMP,NCNT,IOS,NLAY
   CHARACTER(LEN=120) :: FNAME,ONAME
   LOGICAL :: FEXIST

   INTEGER ::IERR,Nsite_tmp
   REAL(DP) :: ODAYS
!   current_time=days2time(time1)
      FNAME = "./"//TRIM(INPUT_DIR)//"/"//trim(casename)//"_el.xy"
!--Make Sure Current Assimilation Data File Exists-----------------------------!
!
   INQUIRE(FILE=TRIM(FNAME),EXIST=FEXIST)
   IF(MSR .AND. .NOT.FEXIST)THEN
     WRITE(IPT,*)'EL OBSERVATION FILE: ',FNAME,' DOES NOT EXIST'
     WRITE(IPT,*)'HALTING.....'
     CALL PSTOP
   END IF
!
!!
!!--Read Number of Current Measurement Stations---------------------------------!
!!
   OPEN(1,FILE=TRIM(FNAME),STATUS='OLD')
   READ(1,*) N_ASSIM_EL

   ALLOCATE(EL0_OBS(N_ASSIM_EL))

!
!--Read X,Y Coordinate of Measurement Stations---------------------------------!
!

   DO I=1,N_ASSIM_EL
     READ(1,*)ITMP,EL0_OBS(I)%X,EL0_OBS(I)%Y !,EL0_OBS(I)%DEPTH,NLAY,CUR_OBS(I)%SITA
     NLAY=1
     EL0_OBS(I)%N_LAYERS = NLAY

     ALLOCATE(EL0_OBS(I)%ODEPTH(NLAY))
     EL0_OBS(I)%chose_flag=0
     DO J=1,NLAY
       READ(1,*)EL0_OBS(I)%ODEPTH(J)
     END DO
   END DO
   CLOSE(1)
!------------------------------------------------------------------------------!
!  Open Current Observation Files for Each Observation Point and Read Data     !
!------------------------------------------------------------------------------!
   ONAME = "./"//TRIM(INPUT_DIR)//"/"//trim(casename)//"_el.dat"
   INQUIRE(FILE=TRIM(ONAME),EXIST=FEXIST)
   IF(MSR .AND. .NOT.FEXIST)THEN
     WRITE(IPT,*)'EL OBSERVATION FILE: ',ONAME,' DOES NOT EXIST'
     WRITE(IPT,*)'HALTING.....'
     CALL PSTOP
   END IF
   OPEN(1,FILE=ONAME,STATUS='old')  ; REWIND(1)

   DO I=1,N_ASSIM_EL
!----Count Number of Data Entries for Observation I----------------------------!

     READ(1,*,IOSTAT=IOS) Nsite_tmp,NCNT     
     IF(IOS<0) then
       WRITE(IPT,*) 'ERROR in read ',trim(casename),'_el.dat at site number:',I
       CALL PSTOP
     END IF
     EL0_OBS(I)%N_TIMES = NCNT
     
     IF(NCNT == 0)THEN
       IF(MSR)WRITE(IPT,*)'NO DATA FOR EL OBSERVATION',I
       CALL PSTOP
     END IF

!----Allocate Arrays to Hold Current (UA,VA) and Time (TIME)-------------------!
     ALLOCATE(EL0_OBS(I)%TIMES(EL0_OBS(I)%N_TIMES))
     ALLOCATE(EL0_OBS(I)%EL(EL0_OBS(I)%N_TIMES , EL0_OBS(I)%N_LAYERS ))
     ALLOCATE(EL0_OBS(I)%obs_days(EL0_OBS(I)%N_TIMES))
!----Read in Current Data for Observation I------------------------------------!
     NLAY = EL0_OBS(I)%N_LAYERS
     DO J=1,EL0_OBS(I)%N_TIMES
       READ(1,*,IOSTAT=IOS) ODAYS,(EL0_OBS(I)%EL(J,K),K=1,NLAY)
       IF(IOS < 0)THEN         
         WRITE(IPT,*)'ERROR in read ',trim(casename),'_el.dat at site:',I,'Time No:',J
         CALL PSTOP
       END IF  
       EL0_OBS(I)%TIMES(J) = DAYS2TIME(ODAYS)
       EL0_OBS(I)%obs_days(J) = ODAYS
       if ( nearequal(current_Time, EL0_OBS(I)%TIMES(J))) then !  .and. ABS(CUR_OBS(I)%UO(J,K)) + ABS(CUR_OBS(I)%VO(J,K)) > .0001) then
           EL0_OBS(I)%chose_flag=1
       end if
     END DO

  END DO
  CLOSE(1)
!---------------------------------------write out obs data for current time
  itmp=0
  IF (MSR) THEN
  ONAME = "rrkf_el.xy"
   OPEN(1,FILE=ONAME,STATUS='replace')  ; REWIND(1)
    write(1,*) sum(EL0_OBS(:)%chose_flag)
   DO I=1,N_ASSIM_EL
     if (EL0_OBS(I)%chose_flag==1) then
     itmp=itmp+1
     write(1,'(i8,3f18.6,i8,f12.6)')ITMP,EL0_OBS(I)%X,EL0_OBS(I)%Y,0.0,EL0_OBS(I)%N_LAYERS
     DO J=1,EL0_OBS(I)%N_LAYERS
       write(1,*)EL0_OBS(I)%ODEPTH(J)
     END DO
     end if
   END DO
   CLOSE(1)
     itmp=0
   ONAME = "rrkf_el.dat"
   OPEN(1,FILE=ONAME,STATUS='replace')  ; REWIND(1)
   DO I=1,N_ASSIM_EL 
     if (EL0_OBS(I)%chose_flag==1) then
         itmp=itmp+1
         write(1,*) itmp,1
         DO J=1,EL0_OBS(I)%N_TIMES
          if (nearequal(EL0_OBS(I)%TIMES(J),current_time)) then
              WRITE(1,'(f18.8,200f12.6)') EL0_OBS(I)%obs_days(j),(EL0_OBS(I)%EL(J,K),K=1,EL0_OBS(I)%N_LAYERS)
          end if
         END DO
     end if
   END DO
  CLOSE(1)
  END IF
   DEALLOCATE(EL0_OBS)
   end subroutine filter_el_obs_data
!-------------------------------------------------------------------------------



   subroutine filter_uv_obs_data(current_time)
   IMPLICIT NONE
   real(dp) :: time1
   TYPE(TIME) :: current_Time
!!-------------------------------------------------
   INTEGER I,J,K,NN,ECNT,ITMP,NCNT,IOS,NLAY
   CHARACTER(LEN=120) :: FNAME,ONAME
   LOGICAL :: FEXIST

   INTEGER ::IERR,Nsite_tmp
   REAL(DP) :: ODAYS
!   current_time=days2time(time1)
      FNAME = "./"//TRIM(INPUT_DIR)//"/"//trim(casename)//"_cur.xy"
!--Make Sure Current Assimilation Data File Exists-----------------------------!
!
   INQUIRE(FILE=TRIM(FNAME),EXIST=FEXIST)
   IF(MSR .AND. .NOT.FEXIST)THEN
     WRITE(IPT,*)'CURRENT OBSERVATION FILE: ',FNAME,' DOES NOT EXIST'
     WRITE(IPT,*)'HALTING.....'
     CALL PSTOP
   END IF
!
!!
!!--Read Number of Current Measurement Stations---------------------------------!
!!
   OPEN(1,FILE=TRIM(FNAME),STATUS='OLD')
   READ(1,*) N_ASSIM_CUR

   ALLOCATE(CUR_OBS(N_ASSIM_CUR))

!
!--Read X,Y Coordinate of Measurement Stations---------------------------------!
!

   DO I=1,N_ASSIM_CUR
     READ(1,*)ITMP,CUR_OBS(I)%X,CUR_OBS(I)%Y,CUR_OBS(I)%DEPTH,NLAY,CUR_OBS(I)%SITA
     CUR_OBS(I)%N_LAYERS = NLAY

     ALLOCATE(CUR_OBS(I)%ODEPTH(NLAY))
     CUR_OBS(I)%chose_flag=0
     DO J=1,NLAY
       READ(1,*)CUR_OBS(I)%ODEPTH(J)
       IF(CUR_OBS(I)%ODEPTH(J) > CUR_OBS(I)%DEPTH)THEN
         IF(MSR)WRITE(IPT,*)'CURRENT OBSERVATION LAYER',J,'OF CURRENT MOORING',I
         IF(MSR)WRITE(IPT,*)'EXCEEDS BATHYMETRIC DEPTH=',CUR_OBS(I)%ODEPTH(J),CUR_OBS(I)%DEPTH
         IF(MSR)WRITE(IPT,*)'HALTING...........'
         CALL PSTOP
       END IF
     END DO
   END DO
   CLOSE(1)
!------------------------------------------------------------------------------!
!  Open Current Observation Files for Each Observation Point and Read Data     !
!------------------------------------------------------------------------------!
   ONAME = "./"//TRIM(INPUT_DIR)//"/"//trim(casename)//"_cur.dat"
   INQUIRE(FILE=TRIM(ONAME),EXIST=FEXIST)
   IF(MSR .AND. .NOT.FEXIST)THEN
     WRITE(IPT,*)'CURRENT OBSERVATION FILE: ',ONAME,' DOES NOT EXIST'
     WRITE(IPT,*)'HALTING.....'
     CALL PSTOP
   END IF
   OPEN(1,FILE=ONAME,STATUS='old')  ; REWIND(1)

   DO I=1,N_ASSIM_CUR
!----Count Number of Data Entries for Observation I----------------------------!

     READ(1,*,IOSTAT=IOS) Nsite_tmp,NCNT     
     IF(IOS<0) then
       WRITE(IPT,*) 'ERROR in read ',trim(casename),'_cur.dat at site number:',I
       CALL PSTOP
     END IF
     CUR_OBS(I)%N_TIMES = NCNT
     
     IF(NCNT == 0)THEN
       IF(MSR)WRITE(IPT,*)'NO DATA FOR CURRENT OBSERVATION',I
       CALL PSTOP
     END IF

!----Allocate Arrays to Hold Current (UA,VA) and Time (TIME)-------------------!
     ALLOCATE(CUR_OBS(I)%TIMES(CUR_OBS(I)%N_TIMES))
     ALLOCATE(CUR_OBS(I)%UO( CUR_OBS(I)%N_TIMES , CUR_OBS(I)%N_LAYERS ))
     ALLOCATE(CUR_OBS(I)%VO( CUR_OBS(I)%N_TIMES , CUR_OBS(I)%N_LAYERS ))
     ALLOCATE(CUR_OBS(I)%obs_days(CUR_OBS(I)%N_TIMES))
!----Read in Current Data for Observation I------------------------------------!
     NLAY = CUR_OBS(I)%N_LAYERS
     DO J=1,CUR_OBS(I)%N_TIMES
       READ(1,*,IOSTAT=IOS) ODAYS,(CUR_OBS(I)%UO(J,K),CUR_OBS(I)%VO(J,K),K=1,NLAY)
       DO K=1,NLAY
          if ( ABS(CUR_OBS(I)%UO(J,K)) + ABS(CUR_OBS(I)%VO(J,K)) < .0001) then
              PRINT *, 'CUR_OBS' ,i, 'VELOCITY', J,K ,'IS ZERO WRONG'
              CALL PSTOP
          end if
       END DO
       IF(IOS < 0)THEN         
         WRITE(IPT,*)'ERROR in read ',trim(casename),'_cur.dat at site:',I,'Time No:',J
         CALL PSTOP
       END IF  
       CUR_OBS(I)%TIMES(J) = DAYS2TIME(ODAYS)
       CUR_OBS(I)%obs_days(J) = ODAYS
       if (nearequal(current_Time, CUR_OBS(I)%TIMES(J))) then !  .and. ABS(CUR_OBS(I)%UO(J,K)) + ABS(CUR_OBS(I)%VO(J,K)) > .0001) then
           CUR_OBS(I)%chose_flag=1
       end if
     END DO

  END DO
  CLOSE(1)
!---------------------------------------write out obs data for current time
  itmp=0
  ONAME = "rrkf_cur.xy"
  IF (MSR) THEN
   OPEN(1,FILE=ONAME,STATUS='replace')  ; REWIND(1)
    write(1,*) sum(CUR_OBS(:)%chose_flag)
   DO I=1,N_ASSIM_CUR
     if (CUR_OBS(I)%chose_flag==1) then
     itmp=itmp+1
     write(1,'(i8,3f18.6,i8,f12.6)')ITMP,CUR_OBS(I)%X,CUR_OBS(I)%Y,CUR_OBS(I)%DEPTH,CUR_OBS(I)%N_LAYERS,CUR_OBS(I)%SITA
     DO J=1,CUR_OBS(I)%N_LAYERS
       write(1,*)CUR_OBS(I)%ODEPTH(J)
     END DO
     end if
   END DO
   CLOSE(1)
     itmp=0
   ONAME = "rrkf_cur.dat"
   OPEN(1,FILE=ONAME,STATUS='replace')  ; REWIND(1)
   DO I=1,N_ASSIM_CUR 
     if (CUR_OBS(I)%chose_flag==1) then
         itmp=itmp+1
         write(1,*) itmp,1
         DO J=1,CUR_OBS(I)%N_TIMES
          if (nearequal(CUR_OBS(I)%TIMES(J),current_time)) then
              WRITE(1,'(f18.8,200f12.6)') CUR_OBS(I)%obs_days(j),(CUR_OBS(I)%UO(J,K),CUR_OBS(I)%VO(J,K),K=1,CUR_OBS(I)%N_LAYERS)
          end if
         END DO
     end if
   END DO
   CLOSE(1)
   END IF
   DEALLOCATE(CUR_OBS)
   end subroutine filter_uv_obs_data
!-------------------------------------------------------------------------------

   subroutine filter_t_obs_data(current_time)
   IMPLICIT NONE
   real(dp) :: time1
   TYPE(TIME) :: current_Time
!!-------------------------------------------------
   INTEGER I,J,K,NN,ECNT,ITMP,NCNT,IOS,NLAY
   CHARACTER(LEN=120) :: FNAME,ONAME
   LOGICAL :: FEXIST

   INTEGER ::IERR,Nsite_tmp
   REAL(DP) :: ODAYS
!   current_time=days2time(time1)
      FNAME = "./"//TRIM(INPUT_DIR)//"/"//trim(casename)//"_t.xy"
!--Make Sure Current Assimilation Data File Exists-----------------------------!
!
   INQUIRE(FILE=TRIM(FNAME),EXIST=FEXIST)
   IF(MSR .AND. .NOT.FEXIST)THEN
     WRITE(IPT,*)'TEMPERATURE OBSERVATION FILE: ',FNAME,' DOES NOT EXIST'
     WRITE(IPT,*)'HALTING.....'
     CALL PSTOP
   END IF
!
!!
!!--Read Number of Current Measurement Stations---------------------------------!
!!
   OPEN(1,FILE=TRIM(FNAME),STATUS='OLD')
   READ(1,*) N_ASSIM_T

   ALLOCATE(T0_OBS(N_ASSIM_T))

!
!--Read X,Y Coordinate of Measurement Stations---------------------------------!
!

   DO I=1,N_ASSIM_T
     READ(1,*)ITMP,T0_OBS(I)%X,T0_OBS(I)%Y,T0_OBS(I)%DEPTH,NLAY,T0_OBS(I)%SITA
     T0_OBS(I)%N_LAYERS = NLAY

     ALLOCATE(T0_OBS(I)%ODEPTH(NLAY))
     T0_OBS(I)%chose_flag=0
     DO J=1,NLAY
       READ(1,*)T0_OBS(I)%ODEPTH(J)
       IF(T0_OBS(I)%ODEPTH(J) > T0_OBS(I)%DEPTH)THEN
         IF(MSR)WRITE(IPT,*)'TEMPERATURE OBSERVATION LAYER',J,'OF TEMPERATURE MOORING',I
         IF(MSR)WRITE(IPT,*)'EXCEEDS BATHYMETRIC DEPTH=',T0_OBS(I)%ODEPTH(J),T0_OBS(I)%DEPTH
         IF(MSR)WRITE(IPT,*)'HALTING...........'
         CALL PSTOP
       END IF
     END DO
   END DO
   CLOSE(1)
!------------------------------------------------------------------------------!
!  Open Current Observation Files for Each Observation Point and Read Data     !
!------------------------------------------------------------------------------!
   ONAME = "./"//TRIM(INPUT_DIR)//"/"//trim(casename)//"_t.dat"
   INQUIRE(FILE=TRIM(ONAME),EXIST=FEXIST)
   IF(MSR .AND. .NOT.FEXIST)THEN
     WRITE(IPT,*)'TEMPERATURE OBSERVATION FILE: ',ONAME,' DOES NOT EXIST'
     WRITE(IPT,*)'HALTING.....'
     CALL PSTOP
   END IF
   OPEN(1,FILE=ONAME,STATUS='old')  ; REWIND(1)

   DO I=1,N_ASSIM_T
!----Count Number of Data Entries for Observation I----------------------------!

     READ(1,*,IOSTAT=IOS) Nsite_tmp,NCNT     
     IF(IOS<0) then
       WRITE(IPT,*) 'ERROR in read ',trim(casename),'_t.dat at site number:',I
       CALL PSTOP
     END IF
     T0_OBS(I)%N_TIMES = NCNT
     
     IF(NCNT == 0)THEN
       IF(MSR)WRITE(IPT,*)'NO DATA FOR CURRENT OBSERVATION',I
       CALL PSTOP
     END IF

!----Allocate Arrays to Hold Current (UA,VA) and Time (TIME)-------------------!
     ALLOCATE(T0_OBS(I)%TIMES(T0_OBS(I)%N_TIMES))
     ALLOCATE(T0_OBS(I)%TEMP( T0_OBS(I)%N_TIMES , T0_OBS(I)%N_LAYERS ))
     ALLOCATE(T0_OBS(I)%obs_days(T0_OBS(I)%N_TIMES))
!----Read in Current Data for Observation I------------------------------------!
     NLAY = T0_OBS(I)%N_LAYERS
     DO J=1,T0_OBS(I)%N_TIMES
       READ(1,*,IOSTAT=IOS) ODAYS,(T0_OBS(I)%TEMP(J,K),K=1,NLAY)
       DO K=1,NLAY
          if ( T0_OBS(I)%TEMP(J,K)< -90.) then
              PRINT *, 'T0_OBS' ,i, 'VELOCITY', J,K ,'IS ZERO WRONG'
              CALL PSTOP
          end if
       END DO
       IF(IOS < 0)THEN         
         WRITE(IPT,*)'ERROR in read ',trim(casename),'_t.dat at site:',I,'Time No:',J
         CALL PSTOP
       END IF  
       T0_OBS(I)%TIMES(J) = DAYS2TIME(ODAYS)
       T0_OBS(I)%obs_days(J) = ODAYS
       if (nearequal(current_Time, T0_OBS(I)%TIMES(J))) then !  .and. ABS(CUR_OBS(I)%UO(J,K)) + ABS(CUR_OBS(I)%VO(J,K)) > .0001) then
           T0_OBS(I)%chose_flag=1
       end if
     END DO

  END DO
  CLOSE(1)
!---------------------------------------write out obs data for current time
  itmp=0
  ONAME = "rrkf_t.xy"
  IF (MSR) THEN
   OPEN(1,FILE=ONAME,STATUS='replace')  ; REWIND(1)
    write(1,*) sum(T0_OBS(:)%chose_flag)
   DO I=1,N_ASSIM_T
     if (T0_OBS(I)%chose_flag==1) then
     itmp=itmp+1
     write(1,'(i8,3f18.6,i8,f12.6)')ITMP,T0_OBS(I)%X,T0_OBS(I)%Y,T0_OBS(I)%DEPTH,T0_OBS(I)%N_LAYERS,T0_OBS(I)%SITA
     DO J=1,T0_OBS(I)%N_LAYERS
       write(1,*)T0_OBS(I)%ODEPTH(J)
     END DO
     end if
   END DO
   CLOSE(1)
     itmp=0
   ONAME = "rrkf_t.dat"
   OPEN(1,FILE=ONAME,STATUS='replace')  ; REWIND(1)
   DO I=1,N_ASSIM_T 
     if (T0_OBS(I)%chose_flag==1) then
         itmp=itmp+1
         write(1,*) itmp,1
         DO J=1,T0_OBS(I)%N_TIMES
          if (nearequal(T0_OBS(I)%TIMES(J),current_time)) then
              WRITE(1,'(f18.8,200f12.6)') T0_OBS(I)%obs_days(j),(T0_OBS(I)%TEMP(J,K),K=1,T0_OBS(I)%N_LAYERS)
          end if
         END DO
     end if
   END DO
   CLOSE(1)
   END IF
   DEALLOCATE(T0_OBS)
   end subroutine filter_t_obs_data
!-------------------------------------------------------------------------------

   subroutine filter_s_obs_data(current_time)
   IMPLICIT NONE
   real(dp) :: time1
   TYPE(TIME) :: current_Time
!!-------------------------------------------------
   INTEGER I,J,K,NN,ECNT,ITMP,NCNT,IOS,NLAY
   CHARACTER(LEN=120) :: FNAME,ONAME
   LOGICAL :: FEXIST

   INTEGER ::IERR,Nsite_tmp
   REAL(DP) :: ODAYS
!   current_time=days2time(time1)
      FNAME = "./"//TRIM(INPUT_DIR)//"/"//trim(casename)//"_s.xy"
!--Make Sure Current Assimilation Data File Exists-----------------------------!
!
   INQUIRE(FILE=TRIM(FNAME),EXIST=FEXIST)
   IF(MSR .AND. .NOT.FEXIST)THEN
     WRITE(IPT,*)'SALINITY OBSERVATION FILE: ',FNAME,' DOES NOT EXIST'
     WRITE(IPT,*)'HALTING.....'
     CALL PSTOP
   END IF
!
!!
!!--Read Number of Current Measurement Stations---------------------------------!
!!
   OPEN(1,FILE=TRIM(FNAME),STATUS='OLD')
   READ(1,*) N_ASSIM_S

   ALLOCATE(S0_OBS(N_ASSIM_S))

!
!--Read X,Y Coordinate of Measurement Stations---------------------------------!
!

   DO I=1,N_ASSIM_S
     READ(1,*)ITMP,S0_OBS(I)%X,S0_OBS(I)%Y,S0_OBS(I)%DEPTH,NLAY,S0_OBS(I)%SITA
     S0_OBS(I)%N_LAYERS = NLAY

     ALLOCATE(S0_OBS(I)%ODEPTH(NLAY))
     S0_OBS(I)%chose_flag=0
     DO J=1,NLAY
       READ(1,*)S0_OBS(I)%ODEPTH(J)
       IF(S0_OBS(I)%ODEPTH(J) > S0_OBS(I)%DEPTH)THEN
         IF(MSR)WRITE(IPT,*)'SALINITY OBSERVATION LAYER',J,'OF SALINITY MOORING',I
         IF(MSR)WRITE(IPT,*)'EXCEEDS BATHYMETRIC DEPTH=',S0_OBS(I)%ODEPTH(J),S0_OBS(I)%DEPTH
         IF(MSR)WRITE(IPT,*)'HALTING...........'
         CALL PSTOP
       END IF
     END DO
   END DO
   CLOSE(1)
!------------------------------------------------------------------------------!
!  Open Current Observation Files for Each Observation Point and Read Data     !
!------------------------------------------------------------------------------!
   ONAME = "./"//TRIM(INPUT_DIR)//"/"//trim(casename)//"_s.dat"
   INQUIRE(FILE=TRIM(ONAME),EXIST=FEXIST)
   IF(MSR .AND. .NOT.FEXIST)THEN
     WRITE(IPT,*)'SALINITY OBSERVATION FILE: ',ONAME,' DOES NOT EXIST'
     WRITE(IPT,*)'HALTING.....'
     CALL PSTOP
   END IF
   OPEN(1,FILE=ONAME,STATUS='old')  ; REWIND(1)

   DO I=1,N_ASSIM_S
!----Count Number of Data Entries for Observation I----------------------------!

     READ(1,*,IOSTAT=IOS) Nsite_tmp,NCNT     
     IF(IOS<0) then
       WRITE(IPT,*) 'ERROR in read ',trim(casename),'_s.dat at site number:',I
       CALL PSTOP
     END IF
     S0_OBS(I)%N_TIMES = NCNT
     
     IF(NCNT == 0)THEN
       IF(MSR)WRITE(IPT,*)'NO DATA FOR SALINITY OBSERVATION',I
       CALL PSTOP
     END IF

!----Allocate Arrays to Hold Current (UA,VA) and Time (TIME)-------------------!
     ALLOCATE(S0_OBS(I)%TIMES(S0_OBS(I)%N_TIMES))
     ALLOCATE(S0_OBS(I)%SAL( S0_OBS(I)%N_TIMES , S0_OBS(I)%N_LAYERS ))
     ALLOCATE(S0_OBS(I)%obs_days(S0_OBS(I)%N_TIMES))
!----Read in Current Data for Observation I------------------------------------!
     NLAY = S0_OBS(I)%N_LAYERS
     DO J=1,S0_OBS(I)%N_TIMES
       READ(1,*,IOSTAT=IOS) ODAYS,(S0_OBS(I)%SAL(J,K),K=1,NLAY)
       DO K=1,NLAY
          if ( S0_OBS(I)%SAL(J,K) <-90) then
              PRINT *, 'SAL_OBS' ,i, 'VELOCITY', J,K ,'IS ZERO WRONG'
              CALL PSTOP
          end if
       END DO
       IF(IOS < 0)THEN         
         WRITE(IPT,*)'ERROR in read ',trim(casename),'_s.dat at site:',I,'Time No:',J
         CALL PSTOP
       END IF  
       S0_OBS(I)%TIMES(J) = DAYS2TIME(ODAYS)
       S0_OBS(I)%obs_days(J) = ODAYS
       if (nearequal(current_Time, S0_OBS(I)%TIMES(J))) then !  .and. ABS(CUR_OBS(I)%UO(J,K)) + ABS(CUR_OBS(I)%VO(J,K)) > .0001) then
           S0_OBS(I)%chose_flag=1
       end if
     END DO

  END DO
  CLOSE(1)
!---------------------------------------write out obs data for current time
  itmp=0
  ONAME = "rrkf_s.xy"
  IF (MSR) THEN
   OPEN(1,FILE=ONAME,STATUS='replace')  ; REWIND(1)
    write(1,*) sum(S0_OBS(:)%chose_flag)
   DO I=1,N_ASSIM_S
     if (S0_OBS(I)%chose_flag==1) then
     itmp=itmp+1
     write(1,'(i8,3f18.6,i8,f12.6)')ITMP,S0_OBS(I)%X,S0_OBS(I)%Y,S0_OBS(I)%DEPTH,S0_OBS(I)%N_LAYERS,S0_OBS(I)%SITA
     DO J=1,S0_OBS(I)%N_LAYERS
       write(1,*)S0_OBS(I)%ODEPTH(J)
     END DO
     end if
   END DO
   CLOSE(1)
     itmp=0
   ONAME = "rrkf_s.dat"
   OPEN(1,FILE=ONAME,STATUS='replace')  ; REWIND(1)
   DO I=1,N_ASSIM_S 
     if (S0_OBS(I)%chose_flag==1) then
         itmp=itmp+1
         write(1,*) itmp,1
         DO J=1,S0_OBS(I)%N_TIMES
          if (NEAREQUAL(S0_OBS(I)%TIMES(J),current_time)) then
              WRITE(1,'(f18.8,200f12.6)') S0_OBS(I)%obs_days(j),(S0_OBS(I)%SAL(J,K),K=1,S0_OBS(I)%N_LAYERS)
          end if
         END DO
     end if
   END DO
   CLOSE(1)
   END IF
   DEALLOCATE(S0_OBS)
   end subroutine filter_s_obs_data
!-------------------------------------------------------------------------------
   subroutine set_RRKF_assim_data
   IMPLICIT NONE
   INTEGER :: I,J,K
   NLOC=0
   IF (EL_ASSIM) THEN
       CALL SET_RRKF_EL_ASSIM_DATA
       NLOC=NLOC+N_ASSIM_EL
   END IF
   IF (UV_ASSIM) THEN
       CALL SET_RRKF_CUR_ASSIM_DATA
      DO I=1,N_ASSIM_CUR
          NLOC=NLOC+2*CUR_OBS(I)%N_LAYERS
      END DO
   END IF
   IF (T_ASSIM) THEN
       CALL SET_RRKF_T_ASSIM_DATA
      DO I=1,N_ASSIM_T
          NLOC=NLOC+T0_OBS(I)%N_LAYERS
      END DO
   END IF
   IF (S_ASSIM) THEN
       print *, 'start to call set_RRKF_s_assim_data'
       CALL SET_RRKF_S_ASSIM_DATA
      DO I=1,N_ASSIM_S
          NLOC=NLOC+S0_OBS(I)%N_LAYERS
      END DO
   END IF
   end subroutine set_RRKF_assim_data
!======================================================================
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%|
   SUBROUTINE SET_RRKF_EL_ASSIM_DATA 
!
!!------------------------------------------------------------------------------!
!!  SET UP ASSIMILATION DATA FOR TEMP/SAL OBSERVATIONS                          |
!!------------------------------------------------------------------------------!
   USE ALL_VARS
#  if defined (MULTIPROCESSOR)   
   USE MOD_PAR  
# endif
   IMPLICIT NONE
   INTEGER I,J,K,NN,ECNT,ITMP,NCNT,IOS,NLAY
   CHARACTER(LEN=120) :: FNAME,ONAME
   CHARACTER(LEN= 2 ) :: NAC   
   INTEGER,  ALLOCATABLE, DIMENSION(:) :: ITEMP
   REAL(SP), ALLOCATABLE, DIMENSION(:) :: FTEMP
   REAL(SP):: X0,Y0,DX,DY,RD,SIGMA_C,ISOBATH_ANGLE,D_ANGLE,ANG_OBS_SIM,DIR_WEIGHT
   REAL(SP), PARAMETER :: ALF = 0.05_SP
   LOGICAL :: FEXIST
   INTEGER :: MAXEL,NBD_CNT
   REAL(SP) :: TS_RADIUS 
   INTEGER :: JMIN,INTMP(1) 
   REAL(SP):: LMIN

   REAL(SP), DIMENSION(1:NGL,1) :: RDLIST
   REAL(SP), DIMENSION(3) :: XTRI,YTRI
   REAL(SP) :: RDLAST
   INTEGER :: LOCIJ(2),MIN_LOC,JJ,IERR,Nsite_tmp
   INTEGER :: ND1,ND2,ND3
   REAL(SP) :: DELTA,COFA,COFB,COFC
   REAL(SP) ::S11,S22,S33,RTMP,RRTMP
   REAL(SP), DIMENSION(KB) :: ZZ_OB
!#  if defined (MULTIPROCESSOR)   
   REAL(SP), ALLOCATABLE :: ZZ_G(:,:)  
!# endif
   REAL(SP) :: X0C,Y0C,HX,HY,H0,TEMP_H
   REAL(DP) :: ODAYS
      
   IF(DBG_SET(DBG_SBR)) write(ipt,*) "START: SET_TS_ASSIM_DATA" 

!------------------------------------------------------------------------------!
!  Read Number of Scalar Observations and Coordinates of Each                  !
!------------------------------------------------------------------------------!
       
   FNAME = "rrkf_el.xy"
!
!--Make Sure elevation Assimilation Data File Exists-----------------------------!
!
   INQUIRE(FILE=TRIM(FNAME),EXIST=FEXIST)
   IF(MSR .AND. .NOT.FEXIST)THEN
     WRITE(IPT,*)'EL OBSERVATION FILE: ',FNAME,' DOES NOT EXIST'
     WRITE(IPT,*)'HALTING.....'
     CALL PSTOP
   END IF
     
!
!--Read Number of T/S Measurement Stations---------------------------------!
!
   OPEN(1,FILE=TRIM(FNAME),STATUS='OLD') ;REWIND(1)
   READ(1,*) N_ASSIM_EL                  !nomber of TS Obs station
   ALLOCATE(EL0_OBS(N_ASSIM_EL))          !Type for TS_OBS

!
!--Read X,Y Coordinate of Measurement Stations---------------------------------!
!
   DO I=1,N_ASSIM_EL  
     READ(1,*)ITMP,EL0_OBS(I)%X,EL0_OBS(I)%Y,EL0_OBS(I)%DEPTH,NLAY,EL0_OBS(I)%SITA
     EL0_OBS(I)%N_LAYERS = NLAY
     ALLOCATE(EL0_OBS(I)%ODEPTH(NLAY))
     DO J=1,NLAY
       READ(1,*)EL0_OBS(I)%ODEPTH(J)
       IF(EL0_OBS(I)%ODEPTH(J) > EL0_OBS(I)%DEPTH)THEN
         IF(MSR)WRITE(IPT,*)'OBSERVATION DEPTH',J,'OF EL OBS',I
         IF(MSR)WRITE(IPT,*)'EXCEEDS BATHYMETRIC DEPTH=',EL0_OBS(I)%ODEPTH(J),EL0_OBS(I)%DEPTH
         IF(MSR)WRITE(IPT,*)'HALTING...........'
         CALL PSTOP
       END IF
     END DO
   END DO

!
!--Shift Coordinates-----------------------------------------------------------!
!
# if !defined (SPHERICAL) 
   EL0_OBS(:)%X = EL0_OBS(:)%X - VXMIN 
   EL0_OBS(:)%Y = EL0_OBS(:)%Y - VYMIN 
# endif
   
!# if defined (SPHERICAL) 
   IF(SERIAL)THEN
     XCG = XC
     YCG = YC
     XG = VX
     YG = VY
   END IF
#  if defined(MULTIPROCESSOR)
   IF(PAR)CALL ACOLLECT(MYID,MSRID,NPROCS,EMAP,XC,XCG)
   IF(PAR)CALL ACOLLECT(MYID,MSRID,NPROCS,EMAP,YC,YCG)
   IF(PAR)CALL ACOLLECT(MYID,MSRID,NPROCS,NMAP,VX,XG)
   IF(PAR)CALL ACOLLECT(MYID,MSRID,NPROCS,NMAP,VY,YG)
   IF(PAR)CALL MPI_BCAST(XCG,NGL+1,MPI_F,0,MPI_FVCOM_GROUP,IERR)
   IF(PAR)CALL MPI_BCAST(YCG,NGL+1,MPI_F,0,MPI_FVCOM_GROUP,IERR)
   IF(PAR)CALL MPI_BCAST(XG,MGL+1,MPI_F,0,MPI_FVCOM_GROUP,IERR)
   IF(PAR)CALL MPI_BCAST(YG,MGL+1,MPI_F,0,MPI_FVCOM_GROUP,IERR)
#  endif
!# endif
!  
!
!--find the cell number (TS_OBS(:)%N_CELL) of Obs station---------------------
!
   RRTMP = 100000.0     !100km
   DO J= 1,N_ASSIM_EL 
     X0 = EL0_OBS(J)%X
     Y0 = EL0_OBS(J)%Y
     DO I=1,NGL
#    if !defined (SPHERICAL)
       Rtmp = SQRT((XCG(I)-X0)*(XCG(I)-X0)+(YCG(I)-Y0)*(YCG(I)-Y0))
#    else
       CALL ARC(X0,Y0,XCG(I),YCG(I),Rtmp)
#    endif
       IF(Rtmp < RRTMP) THEN
	 S11 = (XG(NVG(I,2))-X0)*(YG(NVG(I,3))-Y0)-&
	       (XG(NVG(I,3))-X0)*(YG(NVG(I,2))-Y0)
         S22 = (XG(NVG(I,3))-X0)*(YG(NVG(I,1))-Y0)-&
	       (XG(NVG(I,1))-X0)*(YG(NVG(I,3))-Y0)
         S33 = (XG(NVG(I,1))-X0)*(YG(NVG(I,2))-Y0)-&
	       (XG(NVG(I,2))-X0)*(YG(NVG(I,1))-Y0)
	 IF(S11 <= 0._SP .AND. S22 <= 0._SP .AND. S33 <= 0._SP) THEN
	   EL0_OBS(J)%N_CELL = I
           S11=SQRT((XG(NVG(I,1))-X0)**2+(YG(NVG(I,1))-Y0)**2) ! CACLULATE DISTANCE
           S22=SQRT((XG(NVG(I,2))-X0)**2+(YG(NVG(I,2))-Y0)**2) 
           S33=SQRT((XG(NVG(I,3))-X0)**2+(YG(NVG(I,3))-Y0)**2)
           intmp=minloc([s11,s22,s33])
           EL0_OBS(J)%innode=NVG(I,intmp(1))
	         EXIT
	       ELSE
	         EL0_OBS(J)%N_CELL = 0
	       ENDIF
       ELSE
	 EL0_OBS(J)%N_CELL = -1
       ENDIF
     ENDDO

     IF(EL0_OBS(J)%N_CELL <= 0) THEN
       IF(MSR) WRITE(IPT,*)'ERROR--EL OBS SITE:',J,' OUT OF DOMAN',&
	                    EL0_OBS(J)%N_CELL 
       CALL PSTOP
     ENDIF 	 

   ENDDO             
!--Gather AW0G,AWXG & AWYG use for interp grid to Obs station
   IF(.NOT. ALLOCATED(AW0G)) ALLOCATE(AW0G(0:NGL,3))
   IF(.NOT. ALLOCATED(AWXG)) ALLOCATE(AWXG(0:NGL,3))
   IF(.NOT. ALLOCATED(AWYG)) ALLOCATE(AWYG(0:NGL,3))
   
   IF(SERIAL)THEN
     AW0G = AW0
     AWXG = AWX
     AWYG = AWY
   END IF
#  if defined (MULTIPROCESSOR)
   IF(PAR)THEN     
     CALL ACOLLECT(MYID,MSRID,NPROCS,EMAP,AW0,AW0G)
     CALL ACOLLECT(MYID,MSRID,NPROCS,EMAP,AWX,AWXG)
     CALL ACOLLECT(MYID,MSRID,NPROCS,EMAP,AWY,AWYG)
     CALL MPI_BCAST(AW0G,(NGL+1)*3,MPI_F,0,MPI_FVCOM_GROUP,IERR)
     CALL MPI_BCAST(AWXG,(NGL+1)*3,MPI_F,0,MPI_FVCOM_GROUP,IERR)
     CALL MPI_BCAST(AWYG,(NGL+1)*3,MPI_F,0,MPI_FVCOM_GROUP,IERR)
   END IF
#  endif     
!
!--Close Current Observation Global File---------------------------------------!
!
   CLOSE(1)
       
!------------------------------------------------------------------------------!
!  Open Temp/Sal Observation Files for Each Observation Point and Read Data    !
!------------------------------------------------------------------------------!

!----Make Sure Temperature/Salinity Observation File Exists--------------------!
   ONAME = 'rrkf_el.dat'
   INQUIRE(FILE=TRIM(ONAME),EXIST=FEXIST)
   IF(MSR .AND. .NOT.FEXIST)THEN
     WRITE(IPT,*)'ELEVATION OBSERVATION FILE: ',ONAME,' DOES NOT EXIST'
     WRITE(IPT,*)'HALTING.....'
     CALL PSTOP
   END IF
!----Open Temp/Salinity Observation File for Read------------------------------!
   OPEN(1,FILE=ONAME,STATUS='old')  ; REWIND(1)

   DO I=1,N_ASSIM_EL
     READ(1,*,IOSTAT=IOS) nsite_tmp,NCNT
     IF(IOS<0) then
       WRITE(IPT,*) 'ERROR in read ',trim(casename),'_el.dat at site number:',I
       CALL PSTOP
     ENDIF	
     EL0_OBS(I)%N_TIMES = NCNT
      
     IF(NCNT == 0)THEN
       IF(MSR)WRITE(IPT,*)'NO DATA FOR T OBSERVATION',I
       CALL PSTOP
     END IF

!----Allocate Arrays to Hold EL and Time (TIME)----------!
     ALLOCATE(EL0_OBS(I)%TIMES(EL0_OBS(I)%N_TIMES))
     ALLOCATE(EL0_OBS(I)%EL( EL0_OBS(I)%N_TIMES , EL0_OBS(I)%N_LAYERS ))

!----Read in Current Data for Observation I------------------------------------!
     NLAY = EL0_OBS(I)%N_LAYERS
     DO J=1,EL0_OBS(I)%N_TIMES
       READ(1,*,IOSTAT=IOS)ODAYS,(EL0_OBS(I)%EL(J,K),K=1,NLAY)
       IF(IOS < 0) THEN         ! ios=0 if all goes ok.
         WRITE(IPT,*)'ERROR in read ',trim(casename),'_el.dat at site:',I,  &
                     'Time No:',J
         CALL PSTOP
       ENDIF  
       EL0_OBS(I)%TIMES(J) = DAYS2TIME(ODAYS)
     END DO
  END DO
  CLOSE(1)

!------------------------------------------------------------------------------!
!  Count Number of Points with Bad Data (TEMP = 0. .OR. SAL = 0.)         
!------------------------------------------------------------------------------!
  NBD_CNT = 0
  DO I=1,N_ASSIM_EL
    DO J=1,EL0_OBS(I)%N_TIMES
      DO K=1,EL0_OBS(I)%N_LAYERS
        IF(EL0_OBS(I)%EL(J,K) < -9990.) THEN
          NBD_CNT = NBD_CNT + 1
        END IF
      END DO
    END DO
  END DO
 
!------------------------------------------------------------------------------!
!  Compute Spatial Interpolation Weights for each Mooring Location 
!------------------------------------------------------------------------------!

!------------------------------------------------------------------------------!
!  Compute Sigma Layer Weights for Vertical Interpolation            
!------------------------------------------------------------------------------!
    ALLOCATE(ZZ_G(0:MGL,KB))  
   IF(.NOT. ALLOCATED(HG)) ALLOCATE(HG(0:MGL))
   IF(SERIAL)THEN
     ZZ_G = ZZ
     HG = H
   END IF      
#  if defined (MULTIPROCESSOR)   
   IF(PAR)CALL ACOLLECT(MYID,MSRID,NPROCS,NMAP,ZZ,ZZ_G)
   IF(PAR)CALL ACOLLECT(MYID,MSRID,NPROCS,NMAP,H,HG) 
   IF(PAR)CALL MPI_BCAST(ZZ_G,(MGL+1)*KB,MPI_F,0,MPI_FVCOM_GROUP,IERR)
   IF(PAR)CALL MPI_BCAST(HG,MGL+1,MPI_F,0,MPI_FVCOM_GROUP,IERR)
#  endif

   DO I=1,N_ASSIM_EL
     NLAY = EL0_OBS(I)%N_LAYERS
     ALLOCATE(EL0_OBS(I)%S_INT(NLAY,2))
     ALLOCATE(EL0_OBS(I)%S_WEIGHT(NLAY,2))
     
     X0 = EL0_OBS(I)%X
     Y0 = EL0_OBS(I)%Y
# if defined (SPHERICAL)
     DO NN=1,NGL
       CALL ARC(X0,Y0,XCG(NN),YCG(NN),RDLIST(NN,1))
     END DO
# else
     RDLIST(1:NGL,1) = SQRT((XCG(1:NGL)-X0)**2 + (YCG(1:NGL)-Y0)**2)
# endif
     RDLAST = -1.0_SP
in:  DO WHILE(.TRUE.)
       LOCIJ = MINLOC(RDLIST,RDLIST>RDLAST)
       MIN_LOC = LOCIJ(1)
       IF(MIN_LOC == 0)THEN
         EXIT in
       END IF
       XTRI = XG(NVG(MIN_LOC,1:3))
       YTRI = YG(NVG(MIN_LOC,1:3))
       RDLAST = RDLIST(MIN_LOC,1)
       IF(ISINTRIANGLE1(XTRI,YTRI,X0,Y0))THEN
         JJ = MIN_LOC
	 EXIT IN
       END IF
       RDLAST = RDLIST(MIN_LOC,1)
     END DO IN  	 	      

     ND1=NVG(JJ,1)
     ND2=NVG(JJ,2)
     ND3=NVG(JJ,3)

# if defined (SPHERICAL)
     CALL ARCX(X0,Y0,XCG(JJ),YCG(JJ),X0C)
     CALL ARCX(X0,Y0,XCG(JJ),YCG(JJ),Y0C)
# else
     X0C = X0 - XCG(JJ)
     Y0C = Y0 - YCG(JJ)
# endif
     !----Linear Interpolation of Bathymetry------------------------------------------!
     H0 = AW0G(JJ,1)*HG(ND1)+AW0G(JJ,2)*HG(ND2)+AW0G(JJ,3)*HG(ND3)
     HX = AWXG(JJ,1)*HG(ND1)+AWXG(JJ,2)*HG(ND2)+AWXG(JJ,3)*HG(ND3)
     HY = AWYG(JJ,1)*HG(ND1)+AWYG(JJ,2)*HG(ND2)+AWYG(JJ,3)*HG(ND3)
     TEMP_H = H0 + HX*X0C + HY*Y0C !here is the interpolated observation depth

     DO K=1,KBM1
      !----Linear Interpolation of Bathymetry------------------------------------------!
      H0 = AW0G(JJ,1)*ZZ_G(ND1,K)*HG(ND1)+AW0G(JJ,2)*ZZ_G(ND2,K)*HG(ND2)+AW0G(JJ,3)*ZZ_G(ND3,K)*HG(ND3)
      HX = AWXG(JJ,1)*ZZ_G(ND1,K)*HG(ND1)+AWXG(JJ,2)*ZZ_G(ND2,K)*HG(ND2)+AWXG(JJ,3)*ZZ_G(ND3,K)*HG(ND3)
      HY = AWYG(JJ,1)*ZZ_G(ND1,K)*HG(ND1)+AWYG(JJ,2)*ZZ_G(ND2,K)*HG(ND2)+AWYG(JJ,3)*ZZ_G(ND3,K)*HG(ND3)
      ZZ_OB(K) = (H0 + HX*X0C + HY*Y0C)/TEMP_H
     END DO  

     DO J=1,NLAY
       SIGMA_C = -EL0_OBS(I)%ODEPTH(J)/TEMP_H         !TS_OBS(I)%DEPTH
       DO K=2,KBM1
         IF(ZZ_OB(K) <= SIGMA_C .AND. ZZ_OB(K-1) > SIGMA_C)THEN 
           EL0_OBS(I)%S_INT(J,1) = K-1
           EL0_OBS(I)%S_INT(J,2) = K
           EL0_OBS(I)%S_WEIGHT(J,1) = (SIGMA_C-ZZ_OB(K))/(ZZ_OB(K-1)-ZZ_OB(K))
           EL0_OBS(I)%S_WEIGHT(J,2) = 1.0_SP - EL0_OBS(I)%S_WEIGHT(J,1) 
         END IF  
       END DO
       IF(ZZ_OB(1) <= SIGMA_C)THEN  !!OBSERVATION ABOVE CENTROID OF FIRST SIGMA LAYER
         EL0_OBS(I)%S_INT(J,1) = 1
         EL0_OBS(I)%S_INT(J,2) = 1
         EL0_OBS(I)%S_WEIGHT(J,1) = 1.0_SP
         EL0_OBS(I)%S_WEIGHT(J,2) = 0.0_SP
       END IF
       IF(ZZ_OB(KBM1) > SIGMA_C)THEN !!OBSERVATION BELOW CENTROID OF BOTTOM SIGMA LAYER
         EL0_OBS(I)%S_INT(J,1) = KBM1
         EL0_OBS(I)%S_INT(J,2) = KBM1
         EL0_OBS(I)%S_WEIGHT(J,1) = 1.0_SP
         EL0_OBS(I)%S_WEIGHT(J,2) = 0.0_SP
       END IF

     END DO
   END DO

   DEALLOCATE(ZZ_G)  

!------------------------------------------------------------------------------!
!  Report Number of Interpolation Points, Location and Number of Data 
!------------------------------------------------------------------------------!
   IF(.NOT. MSR)RETURN

   WRITE(IPT,*)
   WRITE(IPT,*)'!            TEMP/SALINITY OBSERVATION DATA           '
   WRITE(IPT,*)" MOORING#   X(KM)      Y(KM)  #INTERP PTS  #DATA TIMES  NEAR_NODE  SITA"
   DO I=1,N_ASSIM_EL
!     MAXEL = MAXLOC(EL0_OBS(I)%X_WEIGHT,DIM=1)
!     WRITE(IPT,'(2X,I5,3X,F8.1,3X,F8.1,3X,I6,5X,I6,5X,I6,5X,F8.1)') &
!     I,EL0_OBS(I)%X/1000.,EL0_OBS(I)%Y/1000., &
!       EL0_OBS(I)%N_INTPTS,EL0_OBS(I)%N_TIMES,EL0_OBS(I)%INTPTS(MAXEL),&
!       EL0_OBS(I)%SITA
   END DO
   WRITE(IPT,*)
   WRITE(IPT,*)'NUMBER OF BAD TS DATA POINTS: ',NBD_CNT
   WRITE(IPT,*)" MOORING #   BEGIN TIME  END TIME"
   DO I=1,N_ASSIM_EL
   WRITE(IPT,*)I,EL0_OBS(I)%TIMES(1)/(24.*3600.),&
       EL0_OBS(I)%TIMES(EL0_OBS(I)%N_TIMES)/(24.*3600.)
   END DO

   RETURN
   END SUBROUTINE SET_RRKF_EL_ASSIM_DATA
!==============================================================================|




!==============================================================================!
   SUBROUTINE SET_RRKF_CUR_ASSIM_DATA 
!==============================================================================!

!------------------------------------------------------------------------------!
!  SET UP ASSIMILATION DATA FOR CURRENT OBSERVATIONS                           |
!------------------------------------------------------------------------------!
   USE ALL_VARS
#  if defined (MULTIPROCESSOR)
   USE MOD_PAR
#  endif
   IMPLICIT NONE
   INTEGER I,J,K,NN,ECNT,ITMP,NCNT,IOS,NLAY
   CHARACTER(LEN=120) :: FNAME,ONAME
   CHARACTER(LEN= 2 ) :: NAC   
   INTEGER,  ALLOCATABLE, DIMENSION(:) :: ITEMP
   REAL(SP), ALLOCATABLE, DIMENSION(:) :: FTEMP
   REAL(SP):: X0,Y0,DX,DY,RD,SIGMA_C,ISOBATH_ANGLE,D_ANGLE,ANG_OBS_SIM,DIR_WEIGHT
   REAL(SP), PARAMETER :: ALF = 0.05_SP
   LOGICAL :: FEXIST
   INTEGER :: MAXEL,NBD_CNT
   REAL(SP) :: CUR_RADIUS 
   INTEGER :: JJ

   REAL(SP), DIMENSION(3) :: XTRI,YTRI 
   REAL(SP) :: RDLAST 
   INTEGER :: LOCIJ(2),MIN_LOC,IERR,Nsite_tmp 
   INTEGER :: ND1,ND2,ND3 
   REAL(SP) :: DELTA,COFA,COFB,COFC 
   REAL(SP) ::S11,S22,S33,RTMP,RRTMP 
   REAL(SP), DIMENSION(1:NGL,1) :: RDLIST 
   REAL(SP), DIMENSION(KB) :: ZZ_OB 
   REAL(SP), ALLOCATABLE :: ZZ_G(:,:)   
   REAL(SP) :: X0C,Y0C,HX,HY,H0,TEMP_H
   REAL(DP) :: ODAYS

   IF (DBG_SET(DBG_SBR)) WRITE(IPT,*) "START: SET_CUR_ASSIM_DATA"

!------------------------------------------------------------------------------!
!  Read Number of Current Observations and Coordinates of Each                 !
!------------------------------------------------------------------------------!
    
   FNAME = "rrkf_cur.xy"
!
!--Make Sure Current Assimilation Data File Exists-----------------------------!
!
   INQUIRE(FILE=TRIM(FNAME),EXIST=FEXIST)
   IF(MSR .AND. .NOT.FEXIST)THEN
     WRITE(IPT,*)'CURRENT OBSERVATION FILE: ',FNAME,' DOES NOT EXIST'
     WRITE(IPT,*)'HALTING.....'
     CALL PSTOP
   END IF
     
!
!--Read Number of Current Measurement Stations---------------------------------!
!
   OPEN(1,FILE=TRIM(FNAME),STATUS='OLD') ;rewind(1)
   READ(1,*) N_ASSIM_CUR

   ALLOCATE(CUR_OBS(N_ASSIM_CUR))

!
!--Read X,Y Coordinate of Measurement Stations---------------------------------!
!

   DO I=1,N_ASSIM_CUR
     READ(1,*)ITMP,CUR_OBS(I)%X,CUR_OBS(I)%Y,CUR_OBS(I)%DEPTH,NLAY,CUR_OBS(I)%SITA
     print *, N_ASSIM_CUR,CUR_OBS(I)%X,CUR_OBS(I)%Y,CUR_OBS(I)%DEPTH,NLAY,CUR_OBS(I)%SITA
     CUR_OBS(I)%N_LAYERS = NLAY

     ALLOCATE(CUR_OBS(I)%ODEPTH(NLAY))
     DO J=1,NLAY
       READ(1,*)CUR_OBS(I)%ODEPTH(J)
       IF(CUR_OBS(I)%ODEPTH(J) > CUR_OBS(I)%DEPTH)THEN
         IF(MSR)WRITE(IPT,*)'CURRENT OBSERVATION LAYER',J,'OF CURRENT MOORING',I
         IF(MSR)WRITE(IPT,*)'EXCEEDS BATHYMETRIC DEPTH=',CUR_OBS(I)%ODEPTH(J),CUR_OBS(I)%DEPTH
         IF(MSR)WRITE(IPT,*)'HALTING...........'
         CALL PSTOP
       END IF
     END DO
     DO J=2,NLAY
       IF(CUR_OBS(I)%ODEPTH(J-1) > CUR_OBS(I)%ODEPTH(J))THEN
         IF(MSR)WRITE(IPT,*)'CURRENT OBSERVATION LAYER',J-1,'OF CURRENT MOORING',I
         IF(MSR)WRITE(IPT,*)'EXCEEDS OBSERVATION LAYER',J,'OF CURRENT MOORING',I,CUR_OBS(I)%ODEPTH(J-1),CUR_OBS(I)%ODEPTH(J)
         IF(MSR)WRITE(IPT,*)'HALTING...........'
         CALL PSTOP
       END IF
     END DO
   END DO
   CUR_MAX_LAYER = MAXVAL(CUR_OBS(1:N_ASSIM_CUR)%N_LAYERS)

    PRINT *, 'FINISH READING RRKF_CUR.XY'
IF(PAR)CALL MPI_BARRIER(MPI_FVCOM_GROUP,IERR)
!
!--Shift Coordinates-----------------------------------------------------------!
!
# if !defined (SPHERICAL) 
   CUR_OBS(:)%X = CUR_OBS(:)%X - VXMIN 
   CUR_OBS(:)%Y = CUR_OBS(:)%Y - VYMIN 
# endif

!# if defined (SPHERICAL) 
   IF(SERIAL)THEN
     XCG = XC
     YCG = YC
     XG = VX
     YG = VY
   END IF
#  if defined(MULTIPROCESSOR)
   IF(PAR)CALL ACOLLECT(MYID,MSRID,NPROCS,EMAP,XC,XCG)
   IF(PAR)CALL ACOLLECT(MYID,MSRID,NPROCS,EMAP,YC,YCG)
   IF(PAR)CALL ACOLLECT(MYID,MSRID,NPROCS,NMAP,VX,XG)
   IF(PAR)CALL ACOLLECT(MYID,MSRID,NPROCS,NMAP,VY,YG)
   IF(PAR)CALL MPI_BCAST(XCG,NGL+1,MPI_F,0,MPI_FVCOM_GROUP,IERR)
   IF(PAR)CALL MPI_BCAST(YCG,NGL+1,MPI_F,0,MPI_FVCOM_GROUP,IERR)
   IF(PAR)CALL MPI_BCAST(XG,MGL+1,MPI_F,0,MPI_FVCOM_GROUP,IERR)
   IF(PAR)CALL MPI_BCAST(YG,MGL+1,MPI_F,0,MPI_FVCOM_GROUP,IERR)
#  endif
!# endif
!  


!
!--find cell number (CUR_OBS(J)%N_CELL) of Obs.
!
   RRTMP = 100000.0     !100km
   DO J= 1,N_ASSIM_CUR 
     X0 = CUR_OBS(J)%X
     Y0 = CUR_OBS(J)%Y
     DO I=1,NGL
#    if !defined (SPHERICAL)
       Rtmp = SQRT((XCG(I)-X0)*(XCG(I)-X0)+(YCG(I)-Y0)*(YCG(I)-Y0))
#    else
       CALL ARC(X0,Y0,XCG(I),YCG(I),Rtmp)
#    endif
       IF(Rtmp < RRTMP) THEN
	 S11 = (XG(NVG(I,2))-X0)*(YG(NVG(I,3))-Y0)-&
	       (XG(NVG(I,3))-X0)*(YG(NVG(I,2))-Y0)
         S22 = (XG(NVG(I,3))-X0)*(YG(NVG(I,1))-Y0)-&
	       (XG(NVG(I,1))-X0)*(YG(NVG(I,3))-Y0)
         S33 = (XG(NVG(I,1))-X0)*(YG(NVG(I,2))-Y0)-&
	       (XG(NVG(I,2))-X0)*(YG(NVG(I,1))-Y0)
	 IF(S11 <= 0._SP .AND. S22 <= 0._SP .AND. S33 <= 0._SP) THEN
	   CUR_OBS(J)%N_CELL = I
	   EXIT
	 ELSE
	   CUR_OBS(J)%N_CELL = 0
	 ENDIF
       ELSE
	 CUR_OBS(J)%N_CELL = -1
       ENDIF
     ENDDO

!     IF(MSR) WRITE(200,*) j,X0,Y0,CUR_OBS(J)%N_CELL
     IF(CUR_OBS(J)%N_CELL <= 0) THEN
       IF(MSR) WRITE(IPT,*)'ERROR--CURRENT OBS SITE:',J,' OUT OF DOMAN',&
	       CUR_OBS(J)%N_CELL 
       CALL PSTOP
     ENDIF 	 
   ENDDO      

!--Gather AW0G,AWXG & AWYG use for interp grid to Obs station
   IF(.NOT. ALLOCATED(AW0G))   ALLOCATE(AW0G(0:NGL,3)) 
   IF(.NOT. ALLOCATED(AWXG))   ALLOCATE(AWXG(0:NGL,3))
   IF(.NOT. ALLOCATED(AWYG))   ALLOCATE(AWYG(0:NGL,3))

   IF(SERIAL)THEN
     AW0G = AW0
     AWXG = AWX
     AWYG = AWY
   END IF
# if defined(MULTIPROCESSOR)
   IF(PAR)CALL ACOLLECT(MYID,MSRID,NPROCS,EMAP,AW0,AW0G)
   IF(PAR)CALL ACOLLECT(MYID,MSRID,NPROCS,EMAP,AWX,AWXG)
   IF(PAR)CALL ACOLLECT(MYID,MSRID,NPROCS,EMAP,AWY,AWYG)
   IF(PAR)CALL MPI_BCAST(AW0G,(NGL+1)*3,MPI_F,0,MPI_FVCOM_GROUP,IERR)
   IF(PAR)CALL MPI_BCAST(AWXG,(NGL+1)*3,MPI_F,0,MPI_FVCOM_GROUP,IERR)
   IF(PAR)CALL MPI_BCAST(AWYG,(NGL+1)*3,MPI_F,0,MPI_FVCOM_GROUP,IERR)
# endif
!
!--Close Current Observation Global File---------------------------------------!
!
   CLOSE(1)   

!------------------------------------------------------------------------------!
!  Open Current Observation Files for Each Observation Point and Read Data     !
!------------------------------------------------------------------------------!
   ONAME = "./"//TRIM(INPUT_DIR)//"/"//trim(casename)//"_cur.dat"
   INQUIRE(FILE=TRIM(ONAME),EXIST=FEXIST)
   IF(MSR .AND. .NOT.FEXIST)THEN
     WRITE(IPT,*)'CURRENT OBSERVATION FILE: ',ONAME,' DOES NOT EXIST'
     WRITE(IPT,*)'HALTING.....'
     CALL PSTOP
   END IF
   OPEN(1,FILE=ONAME,STATUS='old')  ; REWIND(1)

   DO I=1,N_ASSIM_CUR
!----Count Number of Data Entries for Observation I----------------------------!

     READ(1,*,IOSTAT=IOS) Nsite_tmp,NCNT     
     IF(IOS<0) then
       WRITE(IPT,*) 'ERROR in read ',trim(casename),'_cur.dat at site number:',I
       CALL PSTOP
     END IF
     CUR_OBS(I)%N_TIMES = NCNT
     
     IF(NCNT == 0)THEN
       IF(MSR)WRITE(IPT,*)'NO DATA FOR CURRENT OBSERVATION',I
       CALL PSTOP
     END IF

!----Allocate Arrays to Hold Current (UA,VA) and Time (TIME)-------------------!
     ALLOCATE(CUR_OBS(I)%TIMES(CUR_OBS(I)%N_TIMES))
     ALLOCATE(CUR_OBS(I)%UO( CUR_OBS(I)%N_TIMES , CUR_OBS(I)%N_LAYERS ))
     ALLOCATE(CUR_OBS(I)%VO( CUR_OBS(I)%N_TIMES , CUR_OBS(I)%N_LAYERS ))

!----Read in Current Data for Observation I------------------------------------!
     NLAY = CUR_OBS(I)%N_LAYERS
     DO J=1,CUR_OBS(I)%N_TIMES
       READ(1,*,IOSTAT=IOS) ODAYS,(CUR_OBS(I)%UO(J,K),CUR_OBS(I)%VO(J,K),K=1,NLAY)
       IF(IOS < 0)THEN         
         WRITE(IPT,*)'ERROR in read ',trim(casename),'_cur.dat at site:',I,'Time No:',J
         CALL PSTOP
       END IF  
       CUR_OBS(I)%TIMES(J) = DAYS2TIME(ODAYS)
     END DO

!----Convert Time to Seconds---------------------------------------------------!
!----Convert Current Data from cm/s to m/s-------------------------------------!
     CUR_OBS(I)%UO = CUR_OBS(I)%UO * .01_SP
     CUR_OBS(I)%VO = CUR_OBS(I)%VO * .01_SP
  END DO
  CLOSE(1)

!------------------------------------------------------------------------------!
!  Count Number of Points with Bad Data (UO = 0. + V0 = 0.)         
!------------------------------------------------------------------------------!
  NBD_CNT = 0
  DO I=1,N_ASSIM_CUR
    DO J=1,CUR_OBS(I)%N_TIMES
      DO K=1,CUR_OBS(I)%N_LAYERS
        IF(ABS(CUR_OBS(I)%UO(J,K)) + ABS(CUR_OBS(I)%VO(J,K)) < .0001) THEN
          NBD_CNT = NBD_CNT + 1
        END IF
      END DO
    END DO
  END DO
 
!------------------------------------------------------------------------------!
!  Compute Spatial Interpolation Weights for each Mooring Location 
!------------------------------------------------------------------------------!
!------------------------------------------------------------------------------!
!  Compute Sigma Layer Weights for Vertical Interpolation                      
!------------------------------------------------------------------------------!
   ALLOCATE(ZZ_G(0:MGL,KB))
   IF(.NOT. ALLOCATED(HG)) ALLOCATE(HG(0:MGL))
   IF(SERIAL)THEN
     ZZ_G = ZZ
     HG = H
   END IF  
#  if defined (MULTIPROCESSOR)   
   IF(PAR)CALL ACOLLECT(MYID,MSRID,NPROCS,NMAP,ZZ,ZZ_G) 
   IF(PAR)CALL ACOLLECT(MYID,MSRID,NPROCS,NMAP,H,HG) 

   IF(PAR)CALL MPI_BCAST(ZZ_G,(MGL+1)*KB,MPI_F,0,MPI_FVCOM_GROUP,IERR)
   IF(PAR)CALL MPI_BCAST(HG,MGL+1,MPI_F,0,MPI_FVCOM_GROUP,IERR)
#  endif
   DO I=1,N_ASSIM_CUR
     NLAY = CUR_OBS(I)%N_LAYERS
     ALLOCATE(CUR_OBS(I)%S_INT(NLAY,2))
     ALLOCATE(CUR_OBS(I)%S_WEIGHT(NLAY,2))
     X0 = CUR_OBS(I)%X
     Y0 = CUR_OBS(I)%Y
# if defined (SPHERICAL)
     DO NN=1,NGL
       CALL ARC(X0,Y0,XCG(NN),YCG(NN),RDLIST(NN,1))
     END DO
# else
     RDLIST(1:NGL,1) = SQRT((XCG(1:NGL)-X0)**2 + (YCG(1:NGL)-Y0)**2)
# endif
     RDLAST = -1.0_SP
in:  DO WHILE(.TRUE.)
       LOCIJ = MINLOC(RDLIST,RDLIST>RDLAST)
       MIN_LOC = LOCIJ(1)
       IF(MIN_LOC == 0)THEN
         EXIT in
       END IF
       XTRI = XG(NVG(MIN_LOC,1:3))
       YTRI = YG(NVG(MIN_LOC,1:3))
       RDLAST = RDLIST(MIN_LOC,1)
       IF(ISINTRIANGLE1(XTRI,YTRI,X0,Y0))THEN
         JJ = MIN_LOC
	 EXIT IN
       END IF
       RDLAST = RDLIST(MIN_LOC,1)
     END DO IN  	 	      

     ND1=NVG(JJ,1)
     ND2=NVG(JJ,2)
     ND3=NVG(JJ,3)

# if defined (SPHERICAL)
     CALL ARCX(X0,Y0,XCG(JJ),YCG(JJ),X0C)
     CALL ARCX(X0,Y0,XCG(JJ),YCG(JJ),Y0C)
# else
     X0C = X0 - XCG(JJ)
     Y0C = Y0 - YCG(JJ)
# endif
     !----Linear Interpolation of Bathymetry------------------------------------------!
     H0 = AW0G(JJ,1)*HG(ND1)+AW0G(JJ,2)*HG(ND2)+AW0G(JJ,3)*HG(ND3)
     HX = AWXG(JJ,1)*HG(ND1)+AWXG(JJ,2)*HG(ND2)+AWXG(JJ,3)*HG(ND3)
     HY = AWYG(JJ,1)*HG(ND1)+AWYG(JJ,2)*HG(ND2)+AWYG(JJ,3)*HG(ND3)
     TEMP_H = H0 + HX*X0C + HY*Y0C !here is the interpolated observation depth

     DO K=1,KBM1
      !----Linear Interpolation of Bathymetry------------------------------------------!
      H0 = AW0G(JJ,1)*ZZ_G(ND1,K)*HG(ND1)+AW0G(JJ,2)*ZZ_G(ND2,K)*HG(ND2)+AW0G(JJ,3)*ZZ_G(ND3,K)*HG(ND3)
      HX = AWXG(JJ,1)*ZZ_G(ND1,K)*HG(ND1)+AWXG(JJ,2)*ZZ_G(ND2,K)*HG(ND2)+AWXG(JJ,3)*ZZ_G(ND3,K)*HG(ND3)
      HY = AWYG(JJ,1)*ZZ_G(ND1,K)*HG(ND1)+AWYG(JJ,2)*ZZ_G(ND2,K)*HG(ND2)+AWYG(JJ,3)*ZZ_G(ND3,K)*HG(ND3)
      ZZ_OB(K) = (H0 + HX*X0C + HY*Y0C)/TEMP_H
     END DO  

     DO J=1,NLAY
       SIGMA_C = -CUR_OBS(I)%ODEPTH(J)/TEMP_H
       DO K=2,KBM1
         IF(ZZ_OB(K) <= SIGMA_C .AND. ZZ_OB(K-1) > SIGMA_C)THEN 
           CUR_OBS(I)%S_INT(J,1) = K-1
           CUR_OBS(I)%S_INT(J,2) = K
           CUR_OBS(I)%S_WEIGHT(J,1) = (SIGMA_C-ZZ_OB(K))/(ZZ_OB(K-1)-ZZ_OB(K))
           CUR_OBS(I)%S_WEIGHT(J,2) = 1.0_SP - CUR_OBS(I)%S_WEIGHT(J,1) 
         END IF  
       END DO
       IF(ZZ_OB(1) <= SIGMA_C)THEN  !!OBSERVATION ABOVE CENTROID OF FIRST SIGMA LAYER
         CUR_OBS(I)%S_INT(J,1) = 1
         CUR_OBS(I)%S_INT(J,2) = 1
         CUR_OBS(I)%S_WEIGHT(J,1) = 1.0_SP
         CUR_OBS(I)%S_WEIGHT(J,2) = 0.0_SP
       END IF
       IF(ZZ_OB(KBM1) >= SIGMA_C)THEN !!OBSERVATION BELOW CENTROID OF BOTTOM SIGMA LAYER
         CUR_OBS(I)%S_INT(J,1) = KBM1
         CUR_OBS(I)%S_INT(J,2) = KBM1
         CUR_OBS(I)%S_WEIGHT(J,1) = 1.0_SP
         CUR_OBS(I)%S_WEIGHT(J,2) = 0.0_SP
       END IF
       if (msr) print *, "obs",i,"lay", j,SIGMA_C,CUR_OBS(I)%S_INT(J,1),CUR_OBS(I)%S_INT(J,2), CUR_OBS(I)%ODEPTH(J),CUR_OBS(I)%DEPTH,'scoor',zz_ob(1:kbm1)
     END DO
   END DO
!#  if defined (MULTIPROCESSOR)   
   DEALLOCATE(ZZ_G)  
!#  endif
!------------------------------------------------------------------------------!
!  Report Number of Interpolation Points, Location and Number of Data 
!------------------------------------------------------------------------------!
   IF(.NOT. MSR)RETURN

   WRITE(IPT,*)
   WRITE(IPT,*)'!            CURRENT OBSERVATION DATA           '
# if defined (SPHERICAL)
   WRITE(IPT,*)" MOORING#    LON       LAT    #INTERP PTS  #DATA TIMES  NEAR_CELL  SITA"
# else
   WRITE(IPT,*)" MOORING#   X(KM)      Y(KM)  #INTERP PTS  #DATA TIMES  NEAR_CELL  SITA"
# endif
   DO I=1,N_ASSIM_CUR
!     MAXEL = MAXLOC(CUR_OBS(I)%X_WEIGHT,DIM=1)
!     WRITE(IPT,'(2X,I5,3X,F8.1,3X,F8.1,3X,I6,5X,I6,5X,I6,5X,F8.1)') &
!# if defined (SPHERICAL)
!     I,CUR_OBS(I)%X,CUR_OBS(I)%Y, &
!# else
!     I,(CUR_OBS(I)%X+VXMIN)/1000.,(CUR_OBS(I)%Y+VYMIN)/1000., &
!# endif
!       CUR_OBS(I)%N_INTPTS,CUR_OBS(I)%N_TIMES,CUR_OBS(I)%INTPTS(MAXEL),&
!       CUR_OBS(I)%SITA
   END DO
   WRITE(IPT,*)
   WRITE(IPT,*)'NUMBER OF BAD CURRENT DATA POINTS: ',NBD_CNT
   WRITE(IPT,*)" MOORING #   BEGIN TIME  END TIME"
   DO I=1,N_ASSIM_CUR
    WRITE(IPT,*)I, CUR_OBS(I)%N_CELL,CUR_OBS(I)%innode,CUR_OBS(I)%s_weight
!   WRITE(IPT,*)I,CUR_OBS(I)%TIMES(1)/(24.*3600.),&
!       CUR_OBS(I)%TIMES(CUR_OBS(I)%N_TIMES)/(24.*3600.)
   WRITE(IPT,*)I,CUR_OBS(I)%TIMES(1)%MJD,&
       CUR_OBS(I)%TIMES(CUR_OBS(I)%N_TIMES)%MJD
   END DO
   IF (DBG_SET(DBG_SBR)) WRITE(IPT,*) "END: SET_CUR_ASSIM_DATA"

 
   RETURN
   END SUBROUTINE SET_RRKF_CUR_ASSIM_DATA

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%|
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%|
   SUBROUTINE SET_RRKF_T_ASSIM_DATA 
!
!!------------------------------------------------------------------------------!
!!  SET UP ASSIMILATION DATA FOR TEMP/SAL OBSERVATIONS                          |
!!------------------------------------------------------------------------------!
   USE ALL_VARS
#  if defined (MULTIPROCESSOR)   
   USE MOD_PAR  
# endif
   IMPLICIT NONE
   INTEGER I,J,K,NN,ECNT,ITMP,NCNT,IOS,NLAY
   CHARACTER(LEN=120) :: FNAME,ONAME
   CHARACTER(LEN= 2 ) :: NAC   
   INTEGER,  ALLOCATABLE, DIMENSION(:) :: ITEMP
   REAL(SP), ALLOCATABLE, DIMENSION(:) :: FTEMP
   REAL(SP):: X0,Y0,DX,DY,RD,SIGMA_C,ISOBATH_ANGLE,D_ANGLE,ANG_OBS_SIM,DIR_WEIGHT
   REAL(SP), PARAMETER :: ALF = 0.05_SP
   LOGICAL :: FEXIST
   INTEGER :: MAXEL,NBD_CNT
   REAL(SP) :: TS_RADIUS 
   INTEGER :: JMIN,INTMP(1) 
   REAL(SP):: LMIN

   REAL(SP), DIMENSION(1:NGL,1) :: RDLIST
   REAL(SP), DIMENSION(3) :: XTRI,YTRI
   REAL(SP) :: RDLAST
   INTEGER :: LOCIJ(2),MIN_LOC,JJ,IERR,Nsite_tmp
   INTEGER :: ND1,ND2,ND3
   REAL(SP) :: DELTA,COFA,COFB,COFC
   REAL(SP) ::S11,S22,S33,RTMP,RRTMP
   REAL(SP), DIMENSION(KB) :: ZZ_OB
!#  if defined (MULTIPROCESSOR)   
   REAL(SP), ALLOCATABLE :: ZZ_G(:,:)  
!# endif
   REAL(SP) :: X0C,Y0C,HX,HY,H0,TEMP_H
   REAL(DP) :: ODAYS
      
   IF(DBG_SET(DBG_SBR)) write(ipt,*) "START: SET_TS_ASSIM_DATA" 

!------------------------------------------------------------------------------!
!  Read Number of Scalar Observations and Coordinates of Each                  !
!------------------------------------------------------------------------------!
       
   FNAME = "rrkf_t.xy"
!
!--Make Sure Temperature and Salinity Assimilation Data File Exists-----------------------------!
!
   INQUIRE(FILE=TRIM(FNAME),EXIST=FEXIST)
   IF(MSR .AND. .NOT.FEXIST)THEN
     WRITE(IPT,*)'TEMP OBSERVATION FILE: ',FNAME,' DOES NOT EXIST'
     WRITE(IPT,*)'HALTING.....'
     CALL PSTOP
   END IF
     
!
!--Read Number of T/S Measurement Stations---------------------------------!
!
   OPEN(1,FILE=TRIM(FNAME),STATUS='OLD')
   READ(1,*) N_ASSIM_T                  !nomber of TS Obs station
   ALLOCATE(T0_OBS(N_ASSIM_T))          !Type for TS_OBS

!
!--Read X,Y Coordinate of Measurement Stations---------------------------------!
!
   DO I=1,N_ASSIM_T  
     READ(1,*)ITMP,T0_OBS(I)%X,T0_OBS(I)%Y,T0_OBS(I)%DEPTH,NLAY,T0_OBS(I)%SITA
     T0_OBS(I)%N_LAYERS = NLAY
     ALLOCATE(T0_OBS(I)%ODEPTH(NLAY))
     DO J=1,NLAY
       READ(1,*)T0_OBS(I)%ODEPTH(J)
       IF(T0_OBS(I)%ODEPTH(J) > T0_OBS(I)%DEPTH)THEN
         IF(MSR)WRITE(IPT,*)'OBSERVATION DEPTH',J,'OF TEMP OBS',I
         IF(MSR)WRITE(IPT,*)'EXCEEDS BATHYMETRIC DEPTH=',T0_OBS(I)%ODEPTH(J),T0_OBS(I)%DEPTH
         IF(MSR)WRITE(IPT,*)'HALTING...........'
         CALL PSTOP
       END IF
     END DO
   END DO

!
!--Shift Coordinates-----------------------------------------------------------!
!
# if !defined (SPHERICAL) 
   T0_OBS(:)%X = T0_OBS(:)%X - VXMIN 
   T0_OBS(:)%Y = T0_OBS(:)%Y - VYMIN 
# endif
   
!# if defined (SPHERICAL) 
   IF(SERIAL)THEN
     XCG = XC
     YCG = YC
     XG = VX
     YG = VY
   END IF
#  if defined(MULTIPROCESSOR)
   IF(PAR)CALL ACOLLECT(MYID,MSRID,NPROCS,EMAP,XC,XCG)
   IF(PAR)CALL ACOLLECT(MYID,MSRID,NPROCS,EMAP,YC,YCG)
   IF(PAR)CALL ACOLLECT(MYID,MSRID,NPROCS,NMAP,VX,XG)
   IF(PAR)CALL ACOLLECT(MYID,MSRID,NPROCS,NMAP,VY,YG)
   IF(PAR)CALL MPI_BCAST(XCG,NGL+1,MPI_F,0,MPI_FVCOM_GROUP,IERR)
   IF(PAR)CALL MPI_BCAST(YCG,NGL+1,MPI_F,0,MPI_FVCOM_GROUP,IERR)
   IF(PAR)CALL MPI_BCAST(XG,MGL+1,MPI_F,0,MPI_FVCOM_GROUP,IERR)
   IF(PAR)CALL MPI_BCAST(YG,MGL+1,MPI_F,0,MPI_FVCOM_GROUP,IERR)
#  endif
!# endif
!  
!
!--find the cell number (TS_OBS(:)%N_CELL) of Obs station---------------------
!
   RRTMP = 100000.0     !100km
   DO J= 1,N_ASSIM_T 
     X0 = T0_OBS(J)%X
     Y0 = T0_OBS(J)%Y
     DO I=1,NGL
#    if !defined (SPHERICAL)
       Rtmp = SQRT((XCG(I)-X0)*(XCG(I)-X0)+(YCG(I)-Y0)*(YCG(I)-Y0))
#    else
       CALL ARC(X0,Y0,XCG(I),YCG(I),Rtmp)
#    endif
       IF(Rtmp < RRTMP) THEN
	 S11 = (XG(NVG(I,2))-X0)*(YG(NVG(I,3))-Y0)-&
	       (XG(NVG(I,3))-X0)*(YG(NVG(I,2))-Y0)
         S22 = (XG(NVG(I,3))-X0)*(YG(NVG(I,1))-Y0)-&
	       (XG(NVG(I,1))-X0)*(YG(NVG(I,3))-Y0)
         S33 = (XG(NVG(I,1))-X0)*(YG(NVG(I,2))-Y0)-&
	       (XG(NVG(I,2))-X0)*(YG(NVG(I,1))-Y0)
	 IF(S11 <= 0._SP .AND. S22 <= 0._SP .AND. S33 <= 0._SP) THEN
	   T0_OBS(J)%N_CELL = I
           S11=SQRT((XG(NVG(I,1))-X0)**2+(YG(NVG(I,1))-Y0)**2) ! CACLULATE DISTANCE
           S22=SQRT((XG(NVG(I,2))-X0)**2+(YG(NVG(I,2))-Y0)**2)
           S33=SQRT((XG(NVG(I,3))-X0)**2+(YG(NVG(I,3))-Y0)**2)
           intmp=minloc([s11,s22,s33])
           T0_OBS(J)%innode=NVG(I,intmp(1))
	   EXIT
	 ELSE
	   T0_OBS(J)%N_CELL = 0
	 ENDIF
       ELSE
	 T0_OBS(J)%N_CELL = -1
       ENDIF
     ENDDO

     IF(T0_OBS(J)%N_CELL <= 0) THEN
       IF(MSR) WRITE(IPT,*)'ERROR--T/S OBS SITE:',J,' OUT OF DOMAN',&
	                    T0_OBS(J)%N_CELL 
       CALL PSTOP
     ENDIF 	 

   ENDDO             
!--Gather AW0G,AWXG & AWYG use for interp grid to Obs station
   IF(.NOT. ALLOCATED(AW0G)) ALLOCATE(AW0G(0:NGL,3))
   IF(.NOT. ALLOCATED(AWXG)) ALLOCATE(AWXG(0:NGL,3))
   IF(.NOT. ALLOCATED(AWYG)) ALLOCATE(AWYG(0:NGL,3))
   
   IF(SERIAL)THEN
     AW0G = AW0
     AWXG = AWX
     AWYG = AWY
   END IF
#  if defined (MULTIPROCESSOR)
   IF(PAR)THEN     
     CALL ACOLLECT(MYID,MSRID,NPROCS,EMAP,AW0,AW0G)
     CALL ACOLLECT(MYID,MSRID,NPROCS,EMAP,AWX,AWXG)
     CALL ACOLLECT(MYID,MSRID,NPROCS,EMAP,AWY,AWYG)
     CALL MPI_BCAST(AW0G,(NGL+1)*3,MPI_F,0,MPI_FVCOM_GROUP,IERR)
     CALL MPI_BCAST(AWXG,(NGL+1)*3,MPI_F,0,MPI_FVCOM_GROUP,IERR)
     CALL MPI_BCAST(AWYG,(NGL+1)*3,MPI_F,0,MPI_FVCOM_GROUP,IERR)
   END IF
#  endif     
!
!--Close Current Observation Global File---------------------------------------!
!
   CLOSE(1)
       
!------------------------------------------------------------------------------!
!  Open Temp/Sal Observation Files for Each Observation Point and Read Data    !
!------------------------------------------------------------------------------!

!----Make Sure Temperature/Salinity Observation File Exists--------------------!
   ONAME = 'rrkf_t.dat'
   INQUIRE(FILE=TRIM(ONAME),EXIST=FEXIST)
   IF(MSR .AND. .NOT.FEXIST)THEN
     WRITE(IPT,*)'TEMP/SALINITY OBSERVATION FILE: ',ONAME,' DOES NOT EXIST'
     WRITE(IPT,*)'HALTING.....'
     CALL PSTOP
   END IF
!----Open Temp/Salinity Observation File for Read------------------------------!
   OPEN(1,FILE=ONAME,STATUS='old')  ; REWIND(1)

   DO I=1,N_ASSIM_T
     READ(1,*,IOSTAT=IOS) nsite_tmp,NCNT
     IF(IOS<0) then
       WRITE(IPT,*) 'ERROR in read ',trim(casename),'_t.dat at site number:',I
       CALL PSTOP
     ENDIF	
     T0_OBS(I)%N_TIMES = NCNT
      
     IF(NCNT == 0)THEN
       IF(MSR)WRITE(IPT,*)'NO DATA FOR T OBSERVATION',I
       CALL PSTOP
     END IF

!----Allocate Arrays to Hold Temp/Salinity (TEMP/SAL) and Time (TIME)----------!
     ALLOCATE(T0_OBS(I)%TIMES(T0_OBS(I)%N_TIMES))
     ALLOCATE(T0_OBS(I)%TEMP( T0_OBS(I)%N_TIMES , T0_OBS(I)%N_LAYERS ))

!----Read in Current Data for Observation I------------------------------------!
     NLAY = T0_OBS(I)%N_LAYERS
     DO J=1,T0_OBS(I)%N_TIMES
       READ(1,*,IOSTAT=IOS)ODAYS,(T0_OBS(I)%TEMP(J,K),K=1,NLAY)
       IF(IOS < 0) THEN         ! ios=0 if all goes ok.
         WRITE(IPT,*)'ERROR in read ',trim(casename),'_t.dat at site:',I,  &
                     'Time No:',J
         CALL PSTOP
       ENDIF  
       T0_OBS(I)%TIMES(J) = DAYS2TIME(ODAYS)
     END DO
  END DO
  CLOSE(1)

!------------------------------------------------------------------------------!
!  Count Number of Points with Bad Data (TEMP = 0. .OR. SAL = 0.)         
!------------------------------------------------------------------------------!
  NBD_CNT = 0
  DO I=1,N_ASSIM_T
    DO J=1,T0_OBS(I)%N_TIMES
      DO K=1,T0_OBS(I)%N_LAYERS
        IF(T0_OBS(I)%TEMP(J,K) < -90.) THEN
          NBD_CNT = NBD_CNT + 1
        END IF
      END DO
    END DO
  END DO
 
!------------------------------------------------------------------------------!
!  Compute Spatial Interpolation Weights for each Mooring Location 
!------------------------------------------------------------------------------!

!------------------------------------------------------------------------------!
!  Compute Sigma Layer Weights for Vertical Interpolation            
!------------------------------------------------------------------------------!
   ALLOCATE(ZZ_G(0:MGL,KB))  
   IF(.NOT. ALLOCATED(HG)) ALLOCATE(HG(0:MGL))
   IF(SERIAL)THEN
     ZZ_G = ZZ
     HG = H
   END IF      
#  if defined (MULTIPROCESSOR)   
   IF(PAR)CALL ACOLLECT(MYID,MSRID,NPROCS,NMAP,ZZ,ZZ_G)
   IF(PAR)CALL ACOLLECT(MYID,MSRID,NPROCS,NMAP,H,HG) 
   IF(PAR)CALL MPI_BCAST(ZZ_G,(MGL+1)*KB,MPI_F,0,MPI_FVCOM_GROUP,IERR)
   IF(PAR)CALL MPI_BCAST(HG,MGL+1,MPI_F,0,MPI_FVCOM_GROUP,IERR)
#  endif

   DO I=1,N_ASSIM_T
     NLAY = T0_OBS(I)%N_LAYERS
     ALLOCATE(T0_OBS(I)%S_INT(NLAY,2))
     ALLOCATE(T0_OBS(I)%S_WEIGHT(NLAY,2))
     
     X0 = T0_OBS(I)%X
     Y0 = T0_OBS(I)%Y
# if defined (SPHERICAL)
     DO NN=1,NGL
       CALL ARC(X0,Y0,XCG(NN),YCG(NN),RDLIST(NN,1))
     END DO
# else
     RDLIST(1:NGL,1) = SQRT((XCG(1:NGL)-X0)**2 + (YCG(1:NGL)-Y0)**2)
# endif
     RDLAST = -1.0_SP
in:  DO WHILE(.TRUE.)
       LOCIJ = MINLOC(RDLIST,RDLIST>RDLAST)
       MIN_LOC = LOCIJ(1)
       IF(MIN_LOC == 0)THEN
         EXIT in
       END IF
       XTRI = XG(NVG(MIN_LOC,1:3))
       YTRI = YG(NVG(MIN_LOC,1:3))
       RDLAST = RDLIST(MIN_LOC,1)
       IF(ISINTRIANGLE1(XTRI,YTRI,X0,Y0))THEN
         JJ = MIN_LOC
	 EXIT IN
       END IF
       RDLAST = RDLIST(MIN_LOC,1)
     END DO IN  	 	      

     ND1=NVG(JJ,1)
     ND2=NVG(JJ,2)
     ND3=NVG(JJ,3)

# if defined (SPHERICAL)
     CALL ARCX(X0,Y0,XCG(JJ),YCG(JJ),X0C)
     CALL ARCX(X0,Y0,XCG(JJ),YCG(JJ),Y0C)
# else
     X0C = X0 - XCG(JJ)
     Y0C = Y0 - YCG(JJ)
# endif
     !----Linear Interpolation of Bathymetry------------------------------------------!
     H0 = AW0G(JJ,1)*HG(ND1)+AW0G(JJ,2)*HG(ND2)+AW0G(JJ,3)*HG(ND3)
     HX = AWXG(JJ,1)*HG(ND1)+AWXG(JJ,2)*HG(ND2)+AWXG(JJ,3)*HG(ND3)
     HY = AWYG(JJ,1)*HG(ND1)+AWYG(JJ,2)*HG(ND2)+AWYG(JJ,3)*HG(ND3)
     TEMP_H = H0 + HX*X0C + HY*Y0C !here is the interpolated observation depth

     DO K=1,KBM1
      !----Linear Interpolation of Bathymetry------------------------------------------!
      H0 = AW0G(JJ,1)*ZZ_G(ND1,K)*HG(ND1)+AW0G(JJ,2)*ZZ_G(ND2,K)*HG(ND2)+AW0G(JJ,3)*ZZ_G(ND3,K)*HG(ND3)
      HX = AWXG(JJ,1)*ZZ_G(ND1,K)*HG(ND1)+AWXG(JJ,2)*ZZ_G(ND2,K)*HG(ND2)+AWXG(JJ,3)*ZZ_G(ND3,K)*HG(ND3)
      HY = AWYG(JJ,1)*ZZ_G(ND1,K)*HG(ND1)+AWYG(JJ,2)*ZZ_G(ND2,K)*HG(ND2)+AWYG(JJ,3)*ZZ_G(ND3,K)*HG(ND3)
      ZZ_OB(K) = (H0 + HX*X0C + HY*Y0C)/TEMP_H
     END DO  

     DO J=1,NLAY
       SIGMA_C = -T0_OBS(I)%ODEPTH(J)/TEMP_H         !TS_OBS(I)%DEPTH
       DO K=2,KBM1
         IF(ZZ_OB(K) <= SIGMA_C .AND. ZZ_OB(K-1) > SIGMA_C)THEN 
           T0_OBS(I)%S_INT(J,1) = K-1
           T0_OBS(I)%S_INT(J,2) = K
           T0_OBS(I)%S_WEIGHT(J,1) = (SIGMA_C-ZZ_OB(K))/(ZZ_OB(K-1)-ZZ_OB(K))
           T0_OBS(I)%S_WEIGHT(J,2) = 1.0_SP - T0_OBS(I)%S_WEIGHT(J,1) 
         END IF  
       END DO
       IF(ZZ_OB(1) <= SIGMA_C)THEN  !!OBSERVATION ABOVE CENTROID OF FIRST SIGMA LAYER
         T0_OBS(I)%S_INT(J,1) = 1
         T0_OBS(I)%S_INT(J,2) = 1
         T0_OBS(I)%S_WEIGHT(J,1) = 1.0_SP
         T0_OBS(I)%S_WEIGHT(J,2) = 0.0_SP
       END IF
       IF(ZZ_OB(KBM1) > SIGMA_C)THEN !!OBSERVATION BELOW CENTROID OF BOTTOM SIGMA LAYER
         T0_OBS(I)%S_INT(J,1) = KBM1
         T0_OBS(I)%S_INT(J,2) = KBM1
         T0_OBS(I)%S_WEIGHT(J,1) = 1.0_SP
         T0_OBS(I)%S_WEIGHT(J,2) = 0.0_SP
       END IF

     END DO
   END DO

   DEALLOCATE(ZZ_G)  

!------------------------------------------------------------------------------!
!  Report Number of Interpolation Points, Location and Number of Data 
!------------------------------------------------------------------------------!
   IF(.NOT. MSR)RETURN

   WRITE(IPT,*)
   WRITE(IPT,*)'!            TEMP/SALINITY OBSERVATION DATA           '
   WRITE(IPT,*)" MOORING#   X(KM)      Y(KM)  #INTERP PTS  #DATA TIMES  NEAR_NODE  SITA"
   DO I=1,N_ASSIM_T
!     MAXEL = MAXLOC(T0_OBS(I)%X_WEIGHT,DIM=1)
!     WRITE(IPT,'(2X,I5,3X,F8.1,3X,F8.1,3X,I6,5X,I6,5X,I6,5X,F8.1)') &
!     I,T0_OBS(I)%X/1000.,T0_OBS(I)%Y/1000., &
!       T0_OBS(I)%N_INTPTS,T0_OBS(I)%N_TIMES,T0_OBS(I)%INTPTS(MAXEL),&
!       T0_OBS(I)%SITA
   END DO
   WRITE(IPT,*)
   WRITE(IPT,*)'NUMBER OF BAD TS DATA POINTS: ',NBD_CNT
   WRITE(IPT,*)" MOORING #   BEGIN TIME  END TIME"
   DO I=1,N_ASSIM_T
   WRITE(IPT,*)I,T0_OBS(I)%TIMES(1)/(24.*3600.),&
       T0_OBS(I)%TIMES(T0_OBS(I)%N_TIMES)/(24.*3600.)
   END DO

   RETURN
   END SUBROUTINE SET_RRKF_T_ASSIM_DATA
!==============================================================================|

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%|
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%|
   SUBROUTINE SET_RRKF_S_ASSIM_DATA 
!
!!------------------------------------------------------------------------------!
!!  SET UP ASSIMILATION DATA FOR TEMP/SAL OBSERVATIONS                          |
!!------------------------------------------------------------------------------!
   USE ALL_VARS
#  if defined (MULTIPROCESSOR)   
   USE MOD_PAR  
# endif
   IMPLICIT NONE
   INTEGER I,J,K,NN,ECNT,ITMP,NCNT,IOS,NLAY
   CHARACTER(LEN=120) :: FNAME,ONAME
   CHARACTER(LEN= 2 ) :: NAC   
   INTEGER,  ALLOCATABLE, DIMENSION(:) :: ITEMP
   REAL(SP), ALLOCATABLE, DIMENSION(:) :: FTEMP
   REAL(SP):: X0,Y0,DX,DY,RD,SIGMA_C,ISOBATH_ANGLE,D_ANGLE,ANG_OBS_SIM,DIR_WEIGHT
   REAL(SP), PARAMETER :: ALF = 0.05_SP
   LOGICAL :: FEXIST
   INTEGER :: MAXEL,NBD_CNT
   REAL(SP) :: TS_RADIUS 
   INTEGER :: JMIN,INTMP(1) 
   REAL(SP):: LMIN

   REAL(SP), DIMENSION(1:NGL,1) :: RDLIST
   REAL(SP), DIMENSION(3) :: XTRI,YTRI
   REAL(SP) :: RDLAST
   INTEGER :: LOCIJ(2),MIN_LOC,JJ,IERR,Nsite_tmp
   INTEGER :: ND1,ND2,ND3
   REAL(SP) :: DELTA,COFA,COFB,COFC
   REAL(SP) ::S11,S22,S33,RTMP,RRTMP
   REAL(SP), DIMENSION(KB) :: ZZ_OB
!#  if defined (MULTIPROCESSOR)   
   REAL(SP), ALLOCATABLE :: ZZ_G(:,:)  
!# endif
   REAL(SP) :: X0C,Y0C,HX,HY,H0,TEMP_H
   REAL(DP) :: ODAYS
      
   IF(DBG_SET(DBG_SBR)) write(ipt,*) "START: SET_TS_ASSIM_DATA" 

!------------------------------------------------------------------------------!
!  Read Number of Scalar Observations and Coordinates of Each                  !
!------------------------------------------------------------------------------!
       
   FNAME = "rrkf_s.xy"
!
!--Make Sure Temperature and Salinity Assimilation Data File Exists-----------------------------!
!
   INQUIRE(FILE=TRIM(FNAME),EXIST=FEXIST)
   IF(MSR .AND. .NOT.FEXIST)THEN
     WRITE(IPT,*)'Salinity OBSERVATION FILE: ',FNAME,' DOES NOT EXIST'
     WRITE(IPT,*)'HALTING.....'
     CALL PSTOP
   END IF
     
!
!--Read Number of T/S Measurement Stations---------------------------------!
!
   OPEN(1,FILE=TRIM(FNAME),STATUS='OLD')
   READ(1,*) N_ASSIM_S                  !nomber of TS Obs station
   ALLOCATE(S0_OBS(N_ASSIM_S))          !Type for TS_OBS

!
!--Read X,Y Coordinate of Measurement Stations---------------------------------!
!
   DO I=1,N_ASSIM_S  
     READ(1,*)ITMP,S0_OBS(I)%X,S0_OBS(I)%Y,S0_OBS(I)%DEPTH,NLAY,S0_OBS(I)%SITA
     S0_OBS(I)%N_LAYERS = NLAY
     ALLOCATE(S0_OBS(I)%ODEPTH(NLAY))
     DO J=1,NLAY
       READ(1,*)S0_OBS(I)%ODEPTH(J)
       IF(S0_OBS(I)%ODEPTH(J) > S0_OBS(I)%DEPTH)THEN
         IF(MSR)WRITE(IPT,*)'OBSERVATION DEPTH',J,'OF salinity OBS',I
         IF(MSR)WRITE(IPT,*)'EXCEEDS BATHYMETRIC DEPTH=',S0_OBS(I)%ODEPTH(J),S0_OBS(I)%DEPTH
         IF(MSR)WRITE(IPT,*)'HALTING...........'
         CALL PSTOP
       END IF
     END DO
   END DO

!
!--Shift Coordinates-----------------------------------------------------------!
!
# if !defined (SPHERICAL) 
   S0_OBS(:)%X = S0_OBS(:)%X - VXMIN 
   S0_OBS(:)%Y = S0_OBS(:)%Y - VYMIN 
# endif
   
!# if defined (SPHERICAL) 
   IF(SERIAL)THEN
     XCG = XC
     YCG = YC
     XG = VX
     YG = VY
   END IF
#  if defined(MULTIPROCESSOR)
   IF(PAR)CALL ACOLLECT(MYID,MSRID,NPROCS,EMAP,XC,XCG)
   IF(PAR)CALL ACOLLECT(MYID,MSRID,NPROCS,EMAP,YC,YCG)
   IF(PAR)CALL ACOLLECT(MYID,MSRID,NPROCS,NMAP,VX,XG)
   IF(PAR)CALL ACOLLECT(MYID,MSRID,NPROCS,NMAP,VY,YG)
   IF(PAR)CALL MPI_BCAST(XCG,NGL+1,MPI_F,0,MPI_FVCOM_GROUP,IERR)
   IF(PAR)CALL MPI_BCAST(YCG,NGL+1,MPI_F,0,MPI_FVCOM_GROUP,IERR)
   IF(PAR)CALL MPI_BCAST(XG,MGL+1,MPI_F,0,MPI_FVCOM_GROUP,IERR)
   IF(PAR)CALL MPI_BCAST(YG,MGL+1,MPI_F,0,MPI_FVCOM_GROUP,IERR)
#  endif
!# endif
!  
!
!--find the cell number (TS_OBS(:)%N_CELL) of Obs station---------------------
!
   RRTMP = 100000.0     !100km
   DO J= 1,N_ASSIM_S 
     X0 = S0_OBS(J)%X
     Y0 = S0_OBS(J)%Y
     DO I=1,NGL
#    if !defined (SPHERICAL)
       Rtmp = SQRT((XCG(I)-X0)*(XCG(I)-X0)+(YCG(I)-Y0)*(YCG(I)-Y0))
#    else
       CALL ARC(X0,Y0,XCG(I),YCG(I),Rtmp)
#    endif
       IF(Rtmp < RRTMP) THEN
	 S11 = (XG(NVG(I,2))-X0)*(YG(NVG(I,3))-Y0)-&
	       (XG(NVG(I,3))-X0)*(YG(NVG(I,2))-Y0)
         S22 = (XG(NVG(I,3))-X0)*(YG(NVG(I,1))-Y0)-&
	       (XG(NVG(I,1))-X0)*(YG(NVG(I,3))-Y0)
         S33 = (XG(NVG(I,1))-X0)*(YG(NVG(I,2))-Y0)-&
	       (XG(NVG(I,2))-X0)*(YG(NVG(I,1))-Y0)
	 IF(S11 <= 0._SP .AND. S22 <= 0._SP .AND. S33 <= 0._SP) THEN
	   S0_OBS(J)%N_CELL = I
           S11=SQRT((XG(NVG(I,1))-X0)**2+(YG(NVG(I,1))-Y0)**2) ! CACLULATE DISTANCE
           S22=SQRT((XG(NVG(I,2))-X0)**2+(YG(NVG(I,2))-Y0)**2) 
           S33=SQRT((XG(NVG(I,3))-X0)**2+(YG(NVG(I,3))-Y0)**2)
           intmp=minloc([s11,s22,s33])
           S0_OBS(J)%innode=NVG(I,intmp(1))
	   EXIT
	 ELSE
	   S0_OBS(J)%N_CELL = 0
	 ENDIF
       ELSE
	 S0_OBS(J)%N_CELL = -1
       ENDIF
     ENDDO

     IF(S0_OBS(J)%N_CELL <= 0) THEN
       IF(MSR) WRITE(IPT,*)'ERROR--S OBS SITE:',J,' OUT OF DOMAN',&
	                    S0_OBS(J)%N_CELL 
       CALL PSTOP
     ENDIF 	 

   ENDDO             
!--Gather AW0G,AWXG & AWYG use for interp grid to Obs station
   IF(.NOT. ALLOCATED(AW0G)) ALLOCATE(AW0G(0:NGL,3))
   IF(.NOT. ALLOCATED(AWXG)) ALLOCATE(AWXG(0:NGL,3))
   IF(.NOT. ALLOCATED(AWYG)) ALLOCATE(AWYG(0:NGL,3))
   
   IF(SERIAL)THEN
     AW0G = AW0
     AWXG = AWX
     AWYG = AWY
   END IF
#  if defined (MULTIPROCESSOR)
   IF(PAR)THEN     
     CALL ACOLLECT(MYID,MSRID,NPROCS,EMAP,AW0,AW0G)
     CALL ACOLLECT(MYID,MSRID,NPROCS,EMAP,AWX,AWXG)
     CALL ACOLLECT(MYID,MSRID,NPROCS,EMAP,AWY,AWYG)
     CALL MPI_BCAST(AW0G,(NGL+1)*3,MPI_F,0,MPI_FVCOM_GROUP,IERR)
     CALL MPI_BCAST(AWXG,(NGL+1)*3,MPI_F,0,MPI_FVCOM_GROUP,IERR)
     CALL MPI_BCAST(AWYG,(NGL+1)*3,MPI_F,0,MPI_FVCOM_GROUP,IERR)
   END IF
#  endif     
!
!--Close Current Observation Global File---------------------------------------!
!
   CLOSE(1)
       
!------------------------------------------------------------------------------!
!  Open Temp/Sal Observation Files for Each Observation Point and Read Data    !
!------------------------------------------------------------------------------!

!----Make Sure Temperature/Salinity Observation File Exists--------------------!
   ONAME = 'rrkf_s.dat'
   INQUIRE(FILE=TRIM(ONAME),EXIST=FEXIST)
   IF(MSR .AND. .NOT.FEXIST)THEN
     WRITE(IPT,*)'TEMP/SALINITY OBSERVATION FILE: ',ONAME,' DOES NOT EXIST'
     WRITE(IPT,*)'HALTING.....'
     CALL PSTOP
   END IF
!----Open Temp/Salinity Observation File for Read------------------------------!
   OPEN(1,FILE=ONAME,STATUS='old')  ; REWIND(1)

   DO I=1,N_ASSIM_S
     READ(1,*,IOSTAT=IOS) nsite_tmp,NCNT
     IF(IOS<0) then
       WRITE(IPT,*) 'ERROR in read ',trim(casename),'_s.dat at site number:',I
       CALL PSTOP
     ENDIF	
     S0_OBS(I)%N_TIMES = NCNT
      
     IF(NCNT == 0)THEN
       IF(MSR)WRITE(IPT,*)'NO DATA FOR S OBSERVATION',I
       CALL PSTOP
     END IF

!----Allocate Arrays to Hold Temp/Salinity (TEMP/SAL) and Time (TIME)----------!
     ALLOCATE(S0_OBS(I)%TIMES(S0_OBS(I)%N_TIMES))
     ALLOCATE(S0_OBS(I)%SAL( S0_OBS(I)%N_TIMES , S0_OBS(I)%N_LAYERS ))

!----Read in Current Data for Observation I------------------------------------!
     NLAY = S0_OBS(I)%N_LAYERS
     DO J=1,S0_OBS(I)%N_TIMES
       READ(1,*,IOSTAT=IOS)ODAYS,(S0_OBS(I)%SAL(J,K),K=1,NLAY)
       IF(IOS < 0) THEN         ! ios=0 if all goes ok.
         WRITE(IPT,*)'ERROR in read ',trim(casename),'_t.dat at site:',I,  &
                     'Time No:',J
         CALL PSTOP
       ENDIF  
       S0_OBS(I)%TIMES(J) = DAYS2TIME(ODAYS)
     END DO
  END DO
  CLOSE(1)

!------------------------------------------------------------------------------!
!  Count Number of Points with Bad Data (TEMP = 0. .OR. SAL = 0.)         
!------------------------------------------------------------------------------!
  NBD_CNT = 0
  DO I=1,N_ASSIM_S
    DO J=1,S0_OBS(I)%N_TIMES
      DO K=1,S0_OBS(I)%N_LAYERS
        IF(S0_OBS(I)%SAL(J,K) < -90.) THEN
          NBD_CNT = NBD_CNT + 1
        END IF
      END DO
    END DO
  END DO
 
!------------------------------------------------------------------------------!
!  Compute Spatial Interpolation Weights for each Mooring Location 
!------------------------------------------------------------------------------!

!------------------------------------------------------------------------------!
!  Compute Sigma Layer Weights for Vertical Interpolation            
!------------------------------------------------------------------------------!
   ALLOCATE(ZZ_G(0:MGL,KB))  
   IF(.NOT. ALLOCATED(HG)) ALLOCATE(HG(0:MGL))
   IF(SERIAL)THEN
     ZZ_G = ZZ
     HG = H
   END IF      
#  if defined (MULTIPROCESSOR)   
   IF(PAR)CALL ACOLLECT(MYID,MSRID,NPROCS,NMAP,ZZ,ZZ_G)
   IF(PAR)CALL ACOLLECT(MYID,MSRID,NPROCS,NMAP,H,HG) 
   IF(PAR)CALL MPI_BCAST(ZZ_G,(MGL+1)*KB,MPI_F,0,MPI_FVCOM_GROUP,IERR)
   IF(PAR)CALL MPI_BCAST(HG,MGL+1,MPI_F,0,MPI_FVCOM_GROUP,IERR)
#  endif

   DO I=1,N_ASSIM_S
     NLAY = S0_OBS(I)%N_LAYERS
     ALLOCATE(S0_OBS(I)%S_INT(NLAY,2))
     ALLOCATE(S0_OBS(I)%S_WEIGHT(NLAY,2))
     
     X0 = S0_OBS(I)%X
     Y0 = S0_OBS(I)%Y
# if defined (SPHERICAL)
     DO NN=1,NGL
       CALL ARC(X0,Y0,XCG(NN),YCG(NN),RDLIST(NN,1))
     END DO
# else
     RDLIST(1:NGL,1) = SQRT((XCG(1:NGL)-X0)**2 + (YCG(1:NGL)-Y0)**2)
# endif
     RDLAST = -1.0_SP
in:  DO WHILE(.TRUE.)
       LOCIJ = MINLOC(RDLIST,RDLIST>RDLAST)
       MIN_LOC = LOCIJ(1)
       IF(MIN_LOC == 0)THEN
         EXIT in
       END IF
       XTRI = XG(NVG(MIN_LOC,1:3))
       YTRI = YG(NVG(MIN_LOC,1:3))
       RDLAST = RDLIST(MIN_LOC,1)
       IF(ISINTRIANGLE1(XTRI,YTRI,X0,Y0))THEN
         JJ = MIN_LOC
	 EXIT IN
       END IF
       RDLAST = RDLIST(MIN_LOC,1)
     END DO IN  	 	      

     ND1=NVG(JJ,1)
     ND2=NVG(JJ,2)
     ND3=NVG(JJ,3)

# if defined (SPHERICAL)
     CALL ARCX(X0,Y0,XCG(JJ),YCG(JJ),X0C)
     CALL ARCX(X0,Y0,XCG(JJ),YCG(JJ),Y0C)
# else
     X0C = X0 - XCG(JJ)
     Y0C = Y0 - YCG(JJ)
# endif
     !----Linear Interpolation of Bathymetry------------------------------------------!
     H0 = AW0G(JJ,1)*HG(ND1)+AW0G(JJ,2)*HG(ND2)+AW0G(JJ,3)*HG(ND3)
     HX = AWXG(JJ,1)*HG(ND1)+AWXG(JJ,2)*HG(ND2)+AWXG(JJ,3)*HG(ND3)
     HY = AWYG(JJ,1)*HG(ND1)+AWYG(JJ,2)*HG(ND2)+AWYG(JJ,3)*HG(ND3)
     TEMP_H = H0 + HX*X0C + HY*Y0C !here is the interpolated observation depth

     DO K=1,KBM1
      !----Linear Interpolation of Bathymetry------------------------------------------!
      H0 = AW0G(JJ,1)*ZZ_G(ND1,K)*HG(ND1)+AW0G(JJ,2)*ZZ_G(ND2,K)*HG(ND2)+AW0G(JJ,3)*ZZ_G(ND3,K)*HG(ND3)
      HX = AWXG(JJ,1)*ZZ_G(ND1,K)*HG(ND1)+AWXG(JJ,2)*ZZ_G(ND2,K)*HG(ND2)+AWXG(JJ,3)*ZZ_G(ND3,K)*HG(ND3)
      HY = AWYG(JJ,1)*ZZ_G(ND1,K)*HG(ND1)+AWYG(JJ,2)*ZZ_G(ND2,K)*HG(ND2)+AWYG(JJ,3)*ZZ_G(ND3,K)*HG(ND3)
      ZZ_OB(K) = (H0 + HX*X0C + HY*Y0C)/TEMP_H
     END DO  

     DO J=1,NLAY
       SIGMA_C = -S0_OBS(I)%ODEPTH(J)/TEMP_H         !TS_OBS(I)%DEPTH
       DO K=2,KBM1
         IF(ZZ_OB(K) <= SIGMA_C .AND. ZZ_OB(K-1) > SIGMA_C)THEN 
           S0_OBS(I)%S_INT(J,1) = K-1
           S0_OBS(I)%S_INT(J,2) = K
           S0_OBS(I)%S_WEIGHT(J,1) = (SIGMA_C-ZZ_OB(K))/(ZZ_OB(K-1)-ZZ_OB(K))
           S0_OBS(I)%S_WEIGHT(J,2) = 1.0_SP - S0_OBS(I)%S_WEIGHT(J,1) 
         END IF  
       END DO
       IF(ZZ_OB(1) <= SIGMA_C)THEN  !!OBSERVATION ABOVE CENTROID OF FIRST SIGMA LAYER
         S0_OBS(I)%S_INT(J,1) = 1
         S0_OBS(I)%S_INT(J,2) = 1
         S0_OBS(I)%S_WEIGHT(J,1) = 1.0_SP
         S0_OBS(I)%S_WEIGHT(J,2) = 0.0_SP
       END IF
       IF(ZZ_OB(KBM1) > SIGMA_C)THEN !!OBSERVATION BELOW CENTROID OF BOTTOM SIGMA LAYER
         S0_OBS(I)%S_INT(J,1) = KBM1
         S0_OBS(I)%S_INT(J,2) = KBM1
         S0_OBS(I)%S_WEIGHT(J,1) = 1.0_SP
         S0_OBS(I)%S_WEIGHT(J,2) = 0.0_SP
       END IF

     END DO
   END DO

   DEALLOCATE(ZZ_G)  

!------------------------------------------------------------------------------!
!  Report Number of Interpolation Points, Location and Number of Data 
!------------------------------------------------------------------------------!
   IF(.NOT. MSR)RETURN

   WRITE(IPT,*)
   WRITE(IPT,*)'!            TEMP/SALINITY OBSERVATION DATA           '
   WRITE(IPT,*)" MOORING#   X(KM)      Y(KM)  #INTERP PTS  #DATA TIMES  NEAR_NODE  SITA"
   DO I=1,N_ASSIM_S
!     MAXEL = MAXLOC(S0_OBS(I)%X_WEIGHT,DIM=1)
!     WRITE(IPT,'(2X,I5,3X,F8.1,3X,F8.1,3X,I6,5X,I6,5X,I6,5X,F8.1)') &
!     I,S0_OBS(I)%X/1000.,S0_OBS(I)%Y/1000., &
!       S0_OBS(I)%N_INTPTS,S0_OBS(I)%N_TIMES,S0_OBS(I)%INTPTS(MAXEL),&
!       S0_OBS(I)%SITA
   END DO
   WRITE(IPT,*)
   WRITE(IPT,*)'NUMBER OF BAD TS DATA POINTS: ',NBD_CNT
   WRITE(IPT,*)" MOORING #   BEGIN TIME  END TIME"
   DO I=1,N_ASSIM_S
   WRITE(IPT,*)I,S0_OBS(I)%TIMES(1)/(24.*3600.),&
       S0_OBS(I)%TIMES(S0_OBS(I)%N_TIMES)/(24.*3600.)
   END DO

   RETURN
   END SUBROUTINE SET_RRKF_S_ASSIM_DATA
!==============================================================================|

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%|
   LOGICAL FUNCTION NEAREQUAL(TIME1,TIME2)
   IMPLICIT NONE
   TYPE(TIME) :: time1,time2
!!-------------------------------------------------
   INTEGER I,J,K
   NEAREQUAL=.FALSE.
   IF (TIME1%MJD==TIME2%MJD .AND. ABS(TIME1%MuSOD-TIME2%MuSOD)<1000*180) then ! that is 3 minute
    NEAREQUAL=.true.
   end if
   END FUNCTION NEAREQUAL


!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%|


!!=============================================================
!!==============================================================================|
   LOGICAL FUNCTION ISINTRIANGLE1(XT,YT,X0,Y0) 
!!==============================================================================|
!  determine if point (x0,y0) is in triangle defined by nodes (xt(3),yt(3))    |
!  using algorithm used for scene rendering in computer graphics               |
!  algorithm works well unless particle happens to lie in a line parallel      |
!  to the edge of a triangle.                                                  |
!  This can cause problems if you use a regular grid, say for idealized        |
!  modelling and you happen to see particles right on edges or parallel to     |
!  edges.                                                                      |

   USE MOD_PREC
   IMPLICIT NONE
   REAL(SP), INTENT(IN) :: X0,Y0
   REAL(SP), INTENT(IN) :: XT(3),YT(3)
   REAL(SP) :: F1,F2,F3
   REAL(SP) :: X1(2)
   REAL(SP) :: X2(2)
   REAL(SP) :: X3(2)
   REAL(SP) :: P(2)

!------------------------------------------------------------------------------|
   IF (DBG_SET(DBG_SBR)) WRITE(IPT,*) "START: ISINTRIANGLE1" 

   ISINTRIANGLE1 = .FALSE.  

   IF(Y0 < MINVAL(YT) .OR. Y0 > MAXVAL(YT))THEN
     ISINTRIANGLE1 = .FALSE.
     RETURN
   END IF
   IF(X0 < MINVAL(XT) .OR. X0 > MAXVAL(XT))THEN
     ISINTRIANGLE1 = .FALSE.
     RETURN
   END IF

   F1 = (Y0-YT(1))*(XT(2)-XT(1)) - (X0-XT(1))*(YT(2)-YT(1))
   F2 = (Y0-YT(3))*(XT(1)-XT(3)) - (X0-XT(3))*(YT(1)-YT(3))
   F3 = (Y0-YT(2))*(XT(3)-XT(2)) - (X0-XT(2))*(YT(3)-YT(2))
  IF(F1*F3 >= 0.0_SP .AND. F3*F2 >= 0.0_SP) ISINTRIANGLE1 = .TRUE.

   IF (DBG_SET(DBG_SBR)) WRITE(IPT,*) "END: ISINTRIANGLE1" 
   RETURN
   END FUNCTION ISINTRIANGLE1
# endif
!==============================================================================|
END MODULE MOD_RRKF_OBS
